{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import the basic libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import torch\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "from IPython import display\n",
    "import os\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize the environment for running the experiment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Paths for loading/storing results.\n",
    "name = 'results/Gaussian_MINEE_MINE' # filename\n",
    "chkpt_name = name+'.pt'              # checkpoint\n",
    "fig_name = name+'.pdf'               # output figure\n",
    "\n",
    "# use GPU if available\n",
    "if torch.cuda.is_available(): \n",
    "    torch.set_default_tensor_type(torch.cuda.FloatTensor)\n",
    "else:\n",
    "    torch.set_default_tensor_type(torch.FloatTensor)\n",
    "\n",
    "# initialize random seed\n",
    "np.random.seed(0)\n",
    "torch.manual_seed(0);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate data using the gaussian model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from data.gaussian import Gaussian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_size = 1000   # sample size\n",
    "rho = 0.9           # model parameter\n",
    "\n",
    "rep = 1             # number of repeated runs\n",
    "d = 6               # number of dimensions for X (and Y)\n",
    "\n",
    "X = np.zeros((rep,sample_size,d))\n",
    "Y = np.zeros((rep,sample_size,d))\n",
    "g = Gaussian(sample_size=sample_size,rho=rho)\n",
    "for i in range(rep):\n",
    "    for j in range(d):\n",
    "        data = g.data\n",
    "        X[i,:,j] = data[:,0]\n",
    "        Y[i,:,j] = data[:,1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A plot of the first dimension of $Y$ against that of $X$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEWCAYAAABv+EDhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2df7SdZXXnvztcco2BVG9hGlHujY2U1sm6S20GZY1VplJDCtWW1o6ty9hx1Szq2KldmjglKCqhC8gsll214zIzUon1FzPE6iLQoONQxmmwJi6MQXCGWG5ACKDXNBAzide7549z9r37PPd5f53znvO+57zfz1pZ3HPP+77P8wbY+3n23s93i6qCEEJI81hW9QQIIYRUAx0AIYQ0FDoAQghpKHQAhBDSUOgACCGkodABEEJIQ6EDIANHRO4RkT8c0Fh/JCJPisizIvKzOa5/REQuHcTc6oCIfFJEtlc9D1INdACkL7QN6cm24X1SRP5aRM4q+Iw1IqIiMtblHM4EcDOA16vqWar6w26ek/J8FZGXlPlMQgYJHQDpJ7+hqmcBeAWAfwXgmgGP/3MAngPggQGPS8hQQAdA+o6qfh/AXQDWhd+JyDIRuUZEZkTkKRHZJSI/0/763vY/j7V3EhdH7h8XkY+IyOPtPx9p/+4XAHzX3f/V2NxE5K3tsX8oItuC7y4SkX0ickxEnhCRj4rI8vZ3Nrdvtef2b0Xk+SJyh4g8LSI/av/8oqS/FxF5n4h8X0SeEZHvisjrssZtf68i8k4R+b/te68TkbXte46LyG1unpeIyGMicrWI/KC9M3tLypyuEJH722P/g4hMZ82XDDGqyj/8U/ofAI8AuLT98/lorcKva3++B8Aftn9+O4CHAfw8gLMA7AbwqfZ3awAogLGUcT4M4D4A/wLAuQD+wY2Tej+AlwJ4FsBrAIyjFS6ac/P+ZQCvAjDWftaDAN7t7lcAL3GffxbAbwN4LoCzAfw3AH+bMPaFAB4FcJ6b69oC434JwCoA/xLAKQD/o/13+DMAvgPgbe1rL2m/083td3wtgBMALmx//0kA29s/vwLAUwBeCeAMAG9r/3scT5sv/wzvH+4ASD/5WxE5BuBrAP4ewJ9HrnkLgJtV9Xuq+iyAPwPw5gJx/7cA+LCqPqWqTwP4EIC35rz3dwDcoar3quopAO8HMG9fquoBVb1PVedU9REAH0fLgEZR1R+q6u2q+mNVfQbA9SnX/xQtw/pSETlTVR9R1cMFxr1RVY+r6gMADgG4u/13+M9o7bZeHlz/flU9pap/D2APgN+NzOkdAD6uql9X1Z+q6q1oOZdXpc2XDC90AKSf/KaqPk9Vp1T1nap6MnLNeQBm3OcZtFa+P5dzjNj95xW491H7oKonACwkikXkF9phnKMichwtB3ZO0sNE5Lki8vF2SOk4WiGs54nIGeG1qvowgHcD+CCAp0TkcyJyXoFxn3Q/n4x89gn3H7XfzUj6O5oC8J52+OdY23mfj9aqP3G+ZHihAyBV8zhahseYRCtk8SRaoY5u7n8859hPoGXgALQMOFphHONjAB4CcIGqrgJwNQBJed570AqVvLJ9/Wvs0bGLVfUzqvrq9vwVwI1djpvF80Vkpfuc9Hf0KIDr207b/jxXVT+bMV8ypNABkKr5LIA/FZEXt8tE/xzA51V1DsDTaIVkfj7j/mtE5FwROQfABwD8Tc6x/zuAK0Tk1e2k6YfR+f/E2QCOA3hWRH4RwB8F9z8ZzO1stFbfx0RkAsC1SQOLyIUi8qsiMg7g/7Xv+2nOcbvhQyKyXER+BcAVaOUnQv4LgKtE5JXSYqWIXC4iZ2fMlwwpdACkam4B8Cm0wiX/hJZx+WMAUNUfoxVH/9/tkMSrIvdvB7AfwEEA3wbwzfbvMmnHz/89gM+gtRv4EYDH3CXvBfD7AJ5Byzh+PnjEBwHc2p7b7wL4CIAVAH6AVmL671KGHwdwQ/vao2glsa/OOW5RjqL1bo8D+DSAq1T1ofAiVd2PVh7go+3rHwbwBznmS4YUUWVDGEJGFRG5BMDfqGpiOSppLtwBEEJIQ6EDIISQhsIQECGENBTuAAghpKF0pbJYFeecc46uWbOm6mkQQshQceDAgR+o6rnh74fKAaxZswb79++vehqEEDJUiMhM7PcMARFCSEOhAyCEkIZCB0AIIQ2FDoAQQhoKHQAhhDQUOgBCCGkodACEEDIAtuzahy279lU9jQ7oAAghpKEM1UEwQggZNmzVf3BmtuPzjk0XVzYno7IdgIg8R0T+UUS+JSIPiMiHqpoLIYQ0kSp3AKcA/KqqPisiZwL4mojcpar3VTgnQggpFVvp12nlb1TmALSlQ/1s++OZ7T/UpiaE1J46GvNuqDQHICJnADgA4CUA/kpVvx65ZjOAzQAwOTk52AkSQkhJ1NFZ1KIhjIg8D8AXAPyxqh5Kum79+vVKNVBCSFWECd3pqQkA9TTuHhE5oKrrw9/XogxUVY8BuAfAZRVPhRBCGkNlISARORfAT1T1mIisAHApgBurmg8hhGRR54RuN1SZA3gBgFvbeYBlAG5T1TsqnA8hhDSKKquADgJ4eVXjE0JItwz7yt+oRQ6AEEKGiTrq+nQDHQAhhDQUagERQkhO6qzr0w3cARBCGkeREM6ohHticAdACCE5YRkoIYQMKUVCOKMW7olBB0AIIQUZFSdQCy2gvFALiBDSK1t27cPho8exdvWqXIZ8FFb+tdYCIoQQEqefSWiGgAghlTHI1XUY07ffZY09zCv/LOgACCGkhgwiCU0HQAgZON0YN39NN8Zw1Eo4y4AOgBBCasggHBYdACFk4BQxbuFu4cqb9uLEqbnc9yeNnUUTdgp0AISQxjIMRr6fc6MDIIRURh7jFtstdGu4u9lx9MNJ1MXx0AEQQhpHE2Qe8sCTwISQkSc0+CvHW2tfyyVMT00A6HQA/Vz52zxi4/aDpJPA3AEQQhrH2tWrOj43beVvcAdACGkM4aq+qtDPoMflDoAQ0heG2Yg2deVv0AEQQkqjbsnUcD51mVdd5kEHQAhJpGjDFJNZHsSc/MEwANi9dUNfxx1FKnMAInI+gF0AVgOYB7BTVf+iqvkQQnrj8NHjOHFqDgdnZivfCbDMMx9V7gDmALxHVb8pImcDOCAiX1bV71Q4J0II8hlQn0i1lb+XWu4XNq6t/HuRhYiR9ZxRciaVOQBVfQLAE+2fnxGRBwG8EAAdACFDhnXXqotxrLrKZ1ioRRmoiKwBcC+Adap6PPhuM4DNADA5OfnLMzMzA58fIU2lqAEt2m6xV7LGKzL/rENaVR3iKoPaloGKyFkAbgfw7tD4A4Cq7gSwE2idAxjw9AghOYjp9NSBYTDOVVLpDkBEzgRwB4C9qnpz1vU8CEZIPYnlAXpdIYeOJCbTUHSsoo1nun1G3ajdDkBEBMAnADyYx/gTQvpLN+Gb0BAfPrpkE59rnDRjT/pHlSGgfw3grQC+LSL3t393tareWeGcCCEpZK1+165e1XMO4PDR4x1NX2JjJzmOrHnnKQttUpP4KquAvgZAqhqfENLiypv24uTpOcy3o8EHZ2Zx5U17cxnxNK3+EFv5m2E/ODOLjdv3YMXysSXGflnJliHPzqSJVJ4EJoTUn25X0N3Gy80ZmWxzljPKmhfVP+PQARDSUMxY2up75fgYTp6ew7rJicIGsmhnr0NHZrFi+diCfIMd6jLCHUFeQimKWB+Ak6fnsGXXPjoB0AEQQnJQ9GBVXr2emBMC8lf1xE4jp91jOYoyGMZqoBA6AEIaSiip0E8xtZOnO+Ua5rVl9GPG2Fbwh48ez7VS97kFy18ArfcJHQSA2ugV1QE6AELIAlkNU/LWxifp9Vjox/DhGr+CL3qYbFA6RMBoCc3RARAywuSRc/ZiamXLOdvKP8SHeMIV/MbtewC0dgl5Vure2Zw8PRcVh+t3r99hhQ6AEAIAS4ywL8+MlXn6n5NWwyuWj3U4Aavu8Y4mDAPNB+IEg+gxUIRREpqjAyBkBKmDnLONaU7ADHuYiLXPK8eTzwPkyQP4+a8cH0t0HMNssMtmWdUTIIQssmXXvkrE1MxQWhXOvLZ+t3J8DNNTEx3O4uDM7JLQzPTUxMJ1YVjnxKk5rFg+hmWChWvM6Nvz7LoQcxp1EpgzwtDSMMIdACEjSJ4wRVLC1zh89HhiDD/Er+jDE7/GiuVjHUY/vG+ZLA3/AIvVQn5+Se8X7mZY7ZMOdwCE1ABb+Xujlbbq7cdOYcemi7F76wasHB/DyvExrF29quNQmK14p6cmFnYK9t3a1as6wi3+Z3vW7q0bloRk1q5ehWXSuuauay7H9NTEQtjHy0GYE6CkQ7lwB0DICJO3mbu/1lb9eZqupz0rdjgrtnL3p4Dt88nTcwvlon4nEdb6p7WDtGu48k+GDoCQGpC3siTJ4ALIdRI2DTOYsTBMOL4ZYtsJJJE2HwvrxMpQLS+wcryVO1ixfCx3kjpW3grQEcSgAyCkYcROx4YrcU+ecwIx6Ya0cI3fIaQRykuHuxEa996gAyCkRhSVXwbQkXBNknJOk2sOV+JAckK225O6aViOIO2wluU8wlV9jFGq0+83TAITUnOKJnyTNHY8XhcnvN5CLlba6ZO+Nh9LyFrS2nPlTXsXGrqcODW38Dl8J18mak1gijqWoo1hSCfcARAyBITCaOGqNkycWpgmpsrpdXmMrNO6IfZ7X86Z97Ru2IAmRvh+RVb1XPnnp9Km8EVhU3jSJMJ6+lioxNi4fc8Sg+qbqfgTsqH+f6zaxsYywlW+xfyB5ORzUvVQWMFjOw77nCYFnSZD4efkoSOoYVN4Qkg2fmXuD0SFidy01TTQaVRDAwx0SjCnNWv35JVXDp/RTbOXPHr/fj4AMiuUCB0AIbUkybDnDbMUMX5pz4xVDIWJ4NipXmBppU4S6yY7V+6xZHDW3Dw2H5aBZkMHQMiQsEyWGkcfIvKVO7ZzMEllHx4J4/dG2tmDWFI5a5cQzm96amJB9M3PJewd4DE5ijwN65N6EJBk6AAIqRnhQSYz7PNabDXrdw+9SCrHdH0sDOX7++Yp0fSOyZ5tz4r18o0Z8Sx9IjaAzw8dACEl0GuYoWgzljQDCXQmVcPksRnutKRynjEAdISFbEfhu375ew8dmY2GtEKBuNjfhU9YZ512ZmlofiqtAhKRWwBcAeApVV2XdT2rgEhdKdMBhMlMC52Eydkw/GPG0eZh5ZY+xh47gRt7rj0nZrRj0gwW1vFGOtw5+BCVXb9764YluQTvmOy7MBTk/07CKiCu+JdS1yqgTwL4KIBdFc+DkK7otT9saMjtOV4JM0Ys/BI7ALZi+ViHIY1JNafp68Sqi+a1Zej9WKFgm30XOq/YqjzmfMK/v6S8BemNSh2Aqt4rImuqnAMhVROLac9rcjes0OmE5ZyhtIM//JVEeDjr4MwslsnSGn1/PqAXaebYO/mWlFklpezxWw5V7wAyEZHNADYDwOTkZMWzIU0lyciknVDNY5h8SMbCJ2nGOik+DuQvEY2R1o0rVOFMWsnbIbOV42NR6Wi7N4myW1KSbGrvAFR1J4CdQCsHUPF0CCmNcCVvK34z6El17L59o4V+rAeu7SbC5u1Z8gu20rdx/YnhcKUfKwtdJovzPXl6ruOwWuwdQvJKPcR+z5V/99TeARBSJXlj/LFGJEUPIlmyNhZaOXRktiOUE2ugDnT20LWdwpZd+zK1d6zz15U37cUyQUcSNzyNnLVKDw92Fa1wSoKhnvKhAyCkItL62hrmDGx1bvhqobWrV3UIuwGdydLDR493hJZiUs+Hjswu0ROy69PCUlaeGVYaGXlj+7F7Sf+p1AGIyGcBXALgHBF5DMC1qvqJKudEiCctNJHUnCTUounFqFm1TSjg5sfxTiLNyBtZukFAp1JoUmze7xT8PIxYxVE3O4FeK61IMlVXAf1eleMTMkiSVvoWetm4fQ/WTU501MPPa2eCNhbK8bX/G67bAwAd1Tx58PX54Wo/3F34e+z5edpDkvrBf2OE5CAtxj+INoVhpU9SC8e0lb9/lpeB9k4llGoA0kNAnljOI9wFdJMHYIev/kEHQEifiTVlARYNpo+9Hz56HMtkMSnrT8L6Z1hi14y4ddTK0zzdV+v4HUd4+Mw+7966obDAms1n7epVS/r6kvpAB0BIm7wrTFvpJzU8ybo/S8xsXhdF0sxw+lg7gAWp6FhVTt4DWjaO/+ydQJh4TsPu847LN5wvQ5eHzqN82BOYkD6zY9PF2LHpYqwcH8PK8TGsm5zAusnFfrt3XXP5wrW7t27A9NTEQqgkXDWvHB/LlIk4eXpu4RRvFqEzioWPtuzah91bNyzMzXoE21xsx7Jieesswsbte3BwZnbBOZmjpAGvH9wBkMbTbZVJ0mnX2POT9H7mtVMHP9baMXaI6+TpOUxPTURDPXadOQD7p63w/XOS8gW+2iiGF4tL+x2pN3QApFF0k0hMq94pIzGZFRKKVf6E4ZsYdo93COsmJxbuiyV8jROn5hacRax+3z6H1Uahg7FSUa7+6wlDQKTxWIjGwhv2uQzSNPVtRb5i+dhClU8s7h6erDVOnm5JNux9/2IIyUIysfCPOQ2bizWZScJ/x1X9aMIdAGkE3YR5sqp3vDRD2mneNGJ1/iGHjsxG6/NNljk2lq/u8Vh+IO1wV8wphLLSvqlMklIoV/71hw6AjDxJPW1D+mGs8vap9W0fQ8xJJBnnsJIHSD4A5u+3OH+enIC9Q3iqNyt8ReoNHQBpBD5BaavZJNJi/uHnpO/C9oY2h7Dbl5G0YvckOQegM9Gbh1if3zRCQ2/vmaYBROoPHQAZWWJaNEDL+NohpW6SwVnEdhu+I5Yd9gpj7HlO8SaR9z7fzjHJCcQSueHz0yQgylL/JP2HSWDSSExkbcuufR3tEq05uVW+eF17f7AJaJVs2u9C52B18sDiYS27pkzDaPH8LB0e393LnNHurRs6avmNok7IDqoBcfVP7gbqC3cAZCRJUuYEkg1wUp7Aeu/6kkg7ietDI/46oGVUY89M64GbpPmftDuYVywZM2y8bljzGBvH/h7CE8B5iJ1BKEv9kwwOOgDSKEJ5ZSO2OvchG298veHzipihkbbVNrB09xBiOw+Pn2NeQTYv8AZ06vh4Tp6e63AcoZS0vVtRfHWRVQWxEqi+0AGQkSSrHDOmn2NhC1s524q+iCaO4Ru6J82h16SpORxvvC22752P4U/3xiSlfbLaOz9v1KenJpa8V1KDdxr/+kMHQEaCvKdzw7JMr1RpRs4bzlBd04dVwlBNbJWepyNWbIUektb+0RNKToSlm3aNx4eqyjLaNP7DAR0AaQThAbAYtroN7zHjuEw6Nf8Pzsx2dMXyjsIkF2I5AEs+l2Ug00JDSXX6PlxlO5zQaGc1dWdbx+GHDoAMNbHTuv50bnhYySQNisS3fdhk4/Y9WLF8DLu3bkhczdv1YbglbOpimFOxbl4x8paI5lH3NPx8fOUTQCPeFOgASGWUIcyWha3OrflJeBDLh2fC+LaFdsKTsmZkY3OwBu1h4tjG9LmH8JBYFnmcQF7HZge48vYOoEMYTegAyFCS53SuGf6i/XH9vTFMwtknS8PvfVjGG27/syWYrTKpbPa+//IlpalAPD6f17FyhzBa0AGQgVOGMFvSwafwFGrWijipXaFfkcccx4lTcwvhJG+8s5K1sV1B0n1pz8mDl2swh1ZUBoOMNnQAZChJU+FMC6ck1bhbYtbIu1uIVdnkJUtuwTdv70Z0zRt9G+/EqbnUhG5Skrzbpjmk3tABkIGTJLKW556N21uJ0tBYhgYq7eSs59CR2US9/SzCUE9ekqSd/ffA4u4kaQw7aJV2ja/w8U4trESigW8miQ5ARO4E8E5VfaRfg4vIZQD+AsAZAP6rqt7Qr7HIaOANGpC+cs2bNPUxcksET08tTRgPClvt+/LT2Hv4XYHp/MSE7wx/MrfowbO8TpuOY7hI2wF8EsDdInIrgJtU9SdlDiwiZwD4KwC/BuAxAN8QkS+p6nfKHIfUlyI7gVDbJxRh85U7afH7rNV3nnaJvZLlUMyQp638DZtvzPjHwka28k9a6bOBe7NIdACqepuI7AHwAQD7ReRTAObd9zf3OPZFAB5W1e8BgIh8DsAbAdABkCixBuQxYqWVtorOG0u3nUY/Vv9eeC3U7skjO5Ek4RxiuYNw19QLWSt/hpCGi6z/2n4C4ASAcQBnwzmAEnghgEfd58cAvDK8SEQ2A9gMAJOTkyUOT6qmiNHYsenihYSrXwHHDoB5zZuQrBV9aFx70ehPwj8vz6o9D6Hj8DuYE6cWy0z9SWZg6S4sKblORpO0HMBlAG4G8CUAr1DVH5c8dkyAdsn/aqq6E8BOAFi/fn0fNuSk7tjpWjOcoWhbSCzuHWu44pOtSSeE+xECyiKcY+gkpqcmlszXi7nZ770zCQ172XST2CfVk7YD2AbgTar6QJ/GfgzA+e7ziwA83qexSA3pxWiYBo9/VugorDomdioXwMK14UGpcJxBOYFl0qrZt/mETdaBpUqd/ndA567KdP5j8w//rmnAm0laDuBX+jz2NwBcICIvBvB9AG8G8Pt9HpMMIRa2MAO/bnJiicHasmvfkgbnWeGUPBr7MePpQ0x5evnmZV7R4azssJnXEgr7CpgjsJLOpBPFg9Lmp+MYLiprCamqcwDeBWAvgAcB3NbH3QapMTs2Xdy14fD3hnHwPKeA81xnWKLVnIbtHnrBnmmtGdMqlGLs2HRxx07It2cEOpPZ4WG3pOfRiDeHSg+CqeqdAO6scg5keLCdQIgZNet+ZYli3yQltkovunL3xtiLqZ04Ndd1qMiv9m1nEeoGhaWtRloS3XZLSX0OCAF4EpiMIGF+oAixpGuMMG8QGv8kAbjw+aHjiP0cqoj6fr5J+INyjOuTJOgAyNCS1gsA6GxIbqEWH/ePGXufdE1b2Set9v1q3eZlY5pjCnvm+rkW2ZWkJW7LNPZ0IKMLHQAZGmKGKEtGOdTJyUoMh/IPWQ1VwtV8GqH4m6/Jj0k0eGfgw0O+pLOIjhIhIXQApHbkbSoeNi4H0HGflzwwfFgG6DTiyyTZQdh1fkcQXmv9g/3K/ODMbMd15mDmNXll7d/Hxrbnx8JT/TLwPN07+tABkJ4YhFEIDdHG7XsWkqOhNpBnx6aLF64NiSl5pq3203YEaSWWJtIW6xucRKxBTLhDoBEmZUAHQGqDrebN2B6cmV0QJ8uT1PUGPa2/rqdo9U7sejPYvhOZVfWYo7Fy0fAQV8yQm8OwZw2ifj8GD4eNPnQApDDWZtCHJGItGcsyGGYAbTXvq2MALEmuFiGP8feSE142okhf3SyRt5jzA5bmFWiESZnQAZBKSKtcMUM/PTXR8X2SsbY2h+aYej2ZG6u/jxn5cJzws4/923O86NrBmdnMg1nzio7rqtwJkNGDDoDkJlylWlLUG+Akrfm8mv9eCuHQkVYIyKplrEmLl07uR3gkdDTesHsRuiIyEjGJao/X4w8lLQjpF3QAZCCEDVxiTsIbd2BR0TK8x7BVtD8YVYY2j6ltdtMLwPccjp0zyIqrh8bfnA1X4aQfiOrwLDXWr1+v+/fvr3oajSerTDOUJADQcfjKG7kwzBIzuuE9wNLTtWldwHrBnIFPzKaNkyQUZ60mPeHfYVjS6qEDIL0gIgdUdX34+8rE4EgzOHl6Llpbb+JnwKJxy6r0CTtpheJs9qzpqQmsHB/rSN5mHdJKwmSVTUzt8NHjC/NdJq33sPGmpyawe+sGrJucWFAstXmEImuhiJsPrx2cme0Yi8af9AuGgEhhsgySzwdk1eCfODW3ZMcQCqKZobTqG1vxhzsGyxmYIqbPPRRR7UzS8UmL3yclc8Pkb0y2IinuT/E20m+4AyADwVbgMaNmq137ft3kxELIxIyircSBReexe+uGjpW+5QzCVfPho8dz9do1/LXrJicWVvo+Uet3BhbmAhZP/xat3PFjDkq7nxDmAEhfCFe8ZqRNshlYNORJ8XFf5x9T27Sm57HvAOCuay4HUKyyJhzLf7YqJD93e7dYDsLeyxvytMqovBIYhBQlKQfAEBAZCGHTEgAdVT/eIIY193mlHMLr8/S/nZ6aSDx0FY7tdymhomeSkS8KjT8ZJHQApGfSDnWlGUOL4x86Mtuh2dNtwjYkrL03DX+/qvfOxodhQueS1JQlbdwYaZLNNPxk0NABkL6RVtZov++m1j4voayCEQsZmayD4Us/Y2GZtHBN0eQttXZIVdABkK7JIxcc1rn7TlpJ9fq2Su8H/qCWzc9W7BaSWjc5ERVxyyLMbdCwk7pDB0BKJ8kxAPmVN8uSQvBnAYBOATmgM/cQk5aOVQ8V1fNJgnr7pGroAEjX+LJI/zlmEC2k0gthhU7e/r0xrILISzPHxsuzmg8rgyjfQIYFOgBSiFgZY5KRDI1rVjvGGGl6/UWMv2+/aIRJYluJZ0lLh44vvD8v1NsnVUMHQAbCusmJqNaPGfgwNh+jrLCQX7EfnJldqNf3Oj5+vKx8BA05GVYqcQAi8iYAHwTwSwAuUlWe7qo5afHqmOEL4+ShcQUWD3MByavoPE3ZV46PYffWDR2tIo1wdxL+HlisSPIJ4bQ5JO18uoFOg1RJVVIQhwBcCeDeisZvDKbJUzVm/E26YeX42EKZpTkJo2gFkMkxeKmIvOGh8CCXOaJwDl7eISlnQOE2MmxUsgNQ1QcBQKRPtX4kF0VWn2lhjtj9tjOwVbo3zrZCj9XL5zH+tnOwuP7G7Xui+YXdWzcs2ankad8IpO88ijRoT7qGFUCkDtQ+ByAimwFsBoDJycmKZzM81MXAJK3E/arbtHqAuApniNXo2zutm5yIqn1u3L5nSYP1mCaQNXT3Gj5h1y9W9pBRpG8OQES+AmB15KttqvrFvM9R1Z0AdgItMbiSptdo8jiHWEVP2jVJO4MkeefQkM7r0vp7O5TlD495/Encjdv3LPl+XjvlmNPkHcKdgZePiDWDCd81JOvvmIljUgf65gBU9dJ+PZsUJ2/ooyzCRKkRE2pLUta0Od91zeULp4iT5BnCKiMf2w9zIJaQTurGVZawGyF1p1I5aBG5B2BMNKIAAAzYSURBVMB781YBUQ46P2lyzOE1MYMXKmSGksihIqaXdo6pZnr8NUCnA7D4flbYxj/Lv4NVAk1PTSTG3Q3vHNLUPXuBDoTUgVrJQYvIbwH4SwDnAtgjIver6oaM20gBwsNKg9an8avvUBAuXJX7lbuFfrwDiom3JTkEE3VLS9KmvXs3Xbho5MmwUlUV0BcAfKGKsckiadU7thJPS4KGOYBDR2Y7+tr6+L+dBbCSzbAZTFJzeN/43XICwGJsvohxj+0IktRKy4JOgdSZ2lcBkd4IJRB6MUimnVPkeu9E1q5eFa3WCROuaZ3ALHQUVgKlhbH6tUKvS6UVId1CB0AWCA1amGy98qa9HYnYsJLFjLLPCfg8wOGjxzuuATqNpe+ta5hWvw8nxU7x9vKeAPvwkmZCB9AQujFstlq3OnkzvBbisdLLsJmKr+mPrfhDwoohf1DMVvw+XJOVsB1UiSVLOcmwQwdAFoipePq2jSEm4hY2fPFhHMP0evIYS3MgfoXuq4Z8JVI3RpeGm5AWdAAkmsz1XbvS5BnmtVMW2ap4/PfAYrzeh2ySDpsZPnQUNmvJY/wHZdjpQMiwQgdAlhCu4P3hq7SafF96mkcCwojV59vvwwNsPn/g7+1lJ5AGdwlklKEDGFGKCJWFVSwWqvFhFm+kfZLXSzsknf61e4wrb9q7JJZv44bESkNZq09IOdABkAWyKmr8d2a88zZpOXl6LrEtZNIhtbBCKKwI8gJuZRt4lniSJkAHUHOKGp60huzhc9J0cPx3oVKm/TN2j5eA9lhuoNsevknNWrLo1pAPWjuJkCqgAxhxslb1Pulrp3ezEqxhS8U0o+oPc/ndQqwBTOz+WM3+ll37+l7JE/6dceVPRhE6gJrS7co1KYGa9pxY2Wb4PKv0iQmoJY1vhPX7QLqej1GGVENRR5EkokfIKML/ukeAWDllrHTSfg7LML3x9/H40FiaFEQRoxpKNIQni4F48tfI6r7Vzcq8yK6hnzpBhFQNHUBNKbpytdO6dl1sxR1KOOQhdBLhOFnzD8kzfmz3U9ZOIO91TPqSJkAHMMSEhjI0zvZ9Vsw+LVHs77dnHDqymCsoQvjctJV/SNr5gV4T5DTypKnQAdScItUqJ07NdTgBb/h7GT8MKYWlnFmtIYtQ5rN6gU6BNAE6gCEmZpz9qtx+DhOasTi63wWEIR4vEeETxr5pS6y+354VG7MXykiQlz0nQoYROoAhJ+y6lVTnH1LUCPqkbYhJPHsJiCJ9A/z1DM8QMjjoACqiTANnxjnr8FLeOvuk+fnmMjFdIC8CF4afyn7fXp5Jp0JICzqAESEtIZvWLatInXu4qg9bRsbmEDsrkDZHrvwJGRx0AAOm7FBHnueljWESC16+ISlXYKTV9idVGPXDoNNJENIbdAAkkywnE9MQIoTUHzqAAVP2yjjteXkPVGWdts0zfjhmLw4h72EzQkhv0AHUkF6cQ3hvmBjuxth347R6PaRFJ0BI/6nEAYjIDgC/AeA0gMMA/p2qHqtiLlVR1DAWrXG33yWViPYyp7R5+rMBRUk71EYIKR9RzdnRo8xBRV4P4KuqOiciNwKAqr4v677169fr/v37+z6/GIOoTgmNaJHm57HqHl+maVU6/Vi9x8Y3rH9A3nH9obYi9xJCkhGRA6q6Pvx9JTsAVb3bfbwPwO9UMY8sqi5JNGOYpbmfh36rWtq8Nm7fA2BR6rlIY5WiOxZCSG/UIQfwdgCfT/pSRDYD2AwAk5OTg5rTAhaGGMQJ1VDOOVT07PberLJOo4wSVZOE6LZ/byhLQQjpH31zACLyFQCrI19tU9Uvtq/ZBmAOwKeTnqOqOwHsBFohoD5MdQmxhGQvz+nGWWQlawe5Oykylp0W9g1kisKVPyGDoW8OQFUvTfteRN4G4AoAr9MqEhEFsMNSgwpL9DJG3gNgafdmyUXngY1UCKk/VSWBLwNwM4DXqurTee8bdBK4W2niWDI3771lPrvbXUJo8PvxHoSQwVGrJDCAjwIYB/BlEQGA+1T1qormkouqjF4voZ5u59zLyp8QMjxUsgPolirLQLuh28NQecI4dc0BEELqR912ACQD6uMTQvoNdwA1IO0AGOPvhJBe4Q5gyOhFoI0QQvJAB1AD2AyFEFIFdAA1h86AENIv6ABqQJUrf+46CGkuy6qeAOmdLbv2sWafEFIY7gAqpMpST5aZEkLoAIaYLCNOo04ISYMOoELKNtSxfr+DGpsQMnzQAQwxsR4A/uwAwzuEkDToAAZMzBj32os37BzWzU6AENI86ABKpFv56F4JO4fxFDEhJA90AAOiX2EZxvIJId1CB1ACoXG/8qa9Cz1x0+rzrd9wv4w2nQEhJA06gAERW6mXeXiLxp4QUhQ6gBJIM+4xw2wrf1bpEEKqhFIQA2bHpovZMJ0QUgvYEKZCuPInhAyCpIYw3AEQQkhDYQ6gQrjyJ4RUCXcAIwjloQkheaADIISQhlJJCEhErgPwRgDzAJ4C8Aeq+ngVcxklKAJHCClCVTuAHao6raovA3AHgA/0e0CGRQghpJNKdgCqetx9XAlgeGpRawx1gQghRaisCkhErgewCcA/A/g3KddtBrAZACYnJwuPw7AIIYTE6ZsDEJGvAFgd+Wqbqn5RVbcB2CYifwbgXQCujT1HVXcC2Am0DoL1a76jBJ0bISQPlZ8EFpEpAHtUdV3Wtb2cBObKnxDSVGp1ElhELnAf3wDgoSrmQQghTaaqHMANInIhWmWgMwCu6veAXPkTQkgnVVUB/XYV4xJCCFmEJ4EJIaSh0AEQQkhDoQMghJCGQgdACCENhQ6AEEIaCh0AIYQ0lMpPAhdBRJ5G69xA2ZwD4Ad9eO6gGYX3GIV3AEbjPUbhHYDReI9e32FKVc8NfzlUDqBfiMj+2DHpYWMU3mMU3gEYjfcYhXcARuM9+vUODAERQkhDoQMghJCGQgfQYmfVEyiJUXiPUXgHYDTeYxTeARiN9+jLOzAHQAghDYU7AEIIaSh0AIQQ0lDoANqIyHUiclBE7heRu0XkvKrnVBQR2SEiD7Xf4wsi8ryq59QNIvImEXlAROZFZKjK90TkMhH5rog8LCL/ser5dIOI3CIiT4nIoarn0i0icr6I/E8RebD939KfVD2nbhCR54jIP4rIt9rv8aFSn88cQAsRWaWqx9s//wcAL1XVvjeqKRMReT2Ar6rqnIjcCACq+r6Kp1UYEfkltJoFfRzAe1W1uz6gA0ZEzgDwfwD8GoDHAHwDwO+p6ncqnVhBROQ1AJ4FsCtPq9Y6IiIvAPACVf2miJwN4ACA3xzCfxcCYKWqPisiZwL4GoA/UdX7yng+dwBtzPi3WQlg6Dyjqt6tqnPtj/cBeFGV8+kWVX1QVb9b9Ty64CIAD6vq91T1NIDPAXhjxXMqjKreC2C26nn0gqo+oarfbP/8DIAHAbyw2lkVR1s82/54ZvtPabaJDsAhIteLyKMA3gLgA1XPp0feDuCuqifRMF4I4FH3+TEModEZNURkDYCXA/h6tTPpDhE5Q0TuB/AUgC+ramnv0SgHICJfEZFDkT9vBABV3aaq5wP4NIB3VTvbOFnv0L5mG4A5tN6jluR5jyFEIr8bup3kKCEiZwG4HcC7g13+0KCqP1XVl6G1o79IREoLy1XVFL4SVPXSnJd+BsAeANf2cTpdkfUOIvI2AFcAeJ3WOMFT4N/FMPEYgPPd5xcBeLyiuTSedsz8dgCfVtXdVc+nV1T1mIjcA+AyAKUk6Bu1A0hDRC5wH98A4KGq5tItInIZgPcBeIOq/rjq+TSQbwC4QEReLCLLAbwZwJcqnlMjaSdPPwHgQVW9uer5dIuInGvVfCKyAsClKNE2sQqojYjcDuBCtKpPZgBcparfr3ZWxRCRhwGMA/hh+1f3DVslEwCIyG8B+EsA5wI4BuB+Vd1Q7azyISK/DuAjAM4AcIuqXl/xlAojIp8FcAlaEsRPArhWVT9R6aQKIiKvBvC/AHwbrf+nAeBqVb2zulkVR0SmAdyK1n9PywDcpqofLu35dACEENJMGAIihJCGQgdACCENhQ6AEEIaCh0AIYQ0FDoAQghpKHQAhHRJW3Hyn0Rkov35+e3PU1XPjZA80AEQ0iWq+iiAjwG4of2rGwDsVNWZ6mZFSH54DoCQHmjLDRwAcAuAdwB4eVsJlJDa0ygtIELKRlV/IiJbAPwdgNfT+JNhgiEgQnpnI4AnAAxl8xTSXOgACOkBEXkZWh3AXgXgT9udqAgZCugACOmStuLkx9DSmj8CYAeA/1TtrAjJDx0AId3zDgBHVPXL7c//GcAvishrK5wTIblhFRAhhDQU7gAIIaSh0AEQQkhDoQMghJCGQgdACCENhQ6AEEIaCh0AIYQ0FDoAQghpKP8fIPVHspjC7+4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(X[0,:,0],Y[0,:,0],label=\"data\",marker=\"+\",color=\"steelblue\")\n",
    "plt.xlabel('X')\n",
    "plt.ylabel('Y')\n",
    "plt.title('Plot of data samples')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize the MINEE model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from model.minee_mine import MINEE_MINE "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "batch_size = 100       # batch size of data sample\n",
    "ref_batch_factor = 10  # batch size expansion factor for reference sample\n",
    "lr = 5e-5              # learning rate\n",
    "\n",
    "minee_mine_list = []\n",
    "for i in range(rep):\n",
    "    minee_mine_list.append(MINEE_MINE(torch.Tensor(X[i]),torch.Tensor(Y[i]),\n",
    "                            batch_size=batch_size,ref_batch_factor=ref_batch_factor,lr=lr))\n",
    "dXY_list = np.zeros((rep,0))\n",
    "dX_list = np.zeros((rep,0))\n",
    "dY_list = np.zeros((rep,0))\n",
    "mine_dXY_list = np.zeros((rep,0))\n",
    "mi_list = np.zeros((rep,0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load previous results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "load_available = True # set to False to prevent loading previous results\n",
    "if load_available and os.path.exists(chkpt_name):\n",
    "    checkpoint = torch.load(\n",
    "        chkpt_name, map_location='cuda' if torch.cuda.is_available() else 'cpu')\n",
    "    dXY_list = checkpoint['dXY_list']\n",
    "    dX_list = checkpoint['dX_list']\n",
    "    dY_list = checkpoint['dY_list']\n",
    "    mine_dXY_list = checkpoint['mine_dXY_list']\n",
    "    mi_list = checkpoint['mi_list']\n",
    "    minee_mine_state_list = checkpoint['minee_mine_state_list']\n",
    "    for i in range(rep):\n",
    "        minee_mine_list[i].load_state_dict(minee_mine_state_list[i])\n",
    "    print('Previous results loaded.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Continuously train the model. The following can be executed repeatedly and after loading previous results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAEICAYAAABGaK+TAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3gU1frA8e+7u0kgCZ1A6EgVBaSKYr1gudd67b0rImDBgo0Smu1iQ5BmF0EU9ScqdqXY6L0ovabSQwLZzZ7fHzMk2WTTN9ls8n6eJ09mzsyceWd38+bsmXLEGINSSqnQ4wh2AEoppUpGE7hSSoUoTeBKKRWiNIErpVSI0gSulFIhShO4UkqFKE3glYCIzBORe4McQ3sRWSEiR0TkoSKsHyci0+3p5iKSKiLOso+08hKRW0Tkh2DHocqPJvAQISLbRSTdTnSJIvKuiEQXs46WImJExFUGIQ4B5hljahhjxhdnQ2PMTmNMtDEmswziqpT8vZfGmI+MMReV0f6C3khQeWkCDy2XG2OigW5AT2BokOPJqQWwLthB5FRG/6iUqjA0gYcgY8we4FugY+5lIuIQkaEiskNEkkTkAxGpZS9eYP8+aLfkzxSRNiIyX0QOiUiKiMzKb78icoWIrBORg3aLrINd/gvwL2CCXW87P9ueZO/niIj8CNTPsSyrNSkiN4rI0lzbDhaROfZ0hIiME5Gd9jeRySJS3V52vojsFpEnRSQBeNcuHyIi8SKyV0TutffVphj1PWa/lvEicleOuKqLyMv2a31IRH7Lse0ZIvKH/VqtEpHzC3hdG4vIZyKSLCLbcnZBicjpIrJURA7b8b1SwHt5p4j8lmNbIyIDRGST/bqPFpHWIvKnXd8nIhJur1tHRL62YzhgTze1l40Fzsnx/k6wy08WkR9FZL+I/C0i1+fY9yUist7e7x4ReTy/41elYIzRnxD4AbYDF9jTzbBau6Pt+XnAvfb03cBmoBUQDXwOfGgvawkYwJWj3pnAs1j/zKsBZ+ez/3bAUeBCIAyry2QzEJ47hny2/xN4BYgAzgWOANNzxwVE2sva5th2CXCjPf0aMAeoC9QAvgKet5edD3iAF+39VAf+DSQAp9p1f2jvq00x6htlH/MlQBpQx14+0T7uJoAT6G3vtwmwz17fYb9m+4AYP6+LA1gGDAfC7fdtK3BxjtftNns6GjijgPfyTuC3HPPGPraa9vEfB36291ELWA/cYa9bD7jGfo1qAJ8C/5ejLp/3F4gCdgF32e9bNyAFONVeHg+cY0/XAboF+2+oMv4EPQD9KeIbZSXwVOAgsAN4E6huL8v647L/QAfk2K494Lb/yPz90X8ATAWaFrL/YcAnOeYdwB7g/Nwx+Nm2uZ0Io3KUzcBPArfnpwPD7em2WAk9EhCsfyKtc9RzJrDNnj4fyACq5Vj+DnZCtufb2PtqU8T60nO9XknAGfbxpwOn+TneJ7H/aeYo+/5EssxV3gvYmavsaeBde3oBMBKon2sdf+/lneRN4GflmF8GPJlj/mXgtXzesy7AgRzzPu8vcAOwMNc2U4AR9vRO4H6gZrD/dirzj3ahhJb/GmNqG2NaGGMGGGPS/azTGCvBn7ADK3k3zKfOIViJbLHdPXJ3Puv51GuM8WK1wJoUIe7GWMngaK648jMDuMmevhmrJZgGxGAl8mV218RB4Du7/IRkY8yxXPvelWM+53RR6ttnjPHkmE/DagnXx/rGssVP/C2A607Uadd7NtAon3Ub51r3GbLfr3uwvv1sFJElInKZnzoKkphjOt3PfDSAiESKyBS7O+gw1j+O2pL/lUEtgF654r4FiLWXX4P1DWSH3XV2ZjHjVkWgJ3kqn71Yf1wnnGj9JuIn2RpjEoD7AETkbOAnEVlgjNnsp95OJ2ZERLC6cvYUIaZ4oI6IROVI4s2xWoj+/ADUF5EuWIl8sF2egpV0TjXWeQB/ctcZDzTNMd8sx3RR6stPCnAMaA2syrVsF1YL/L4i1LMLq8Xf1t9CY8wm4CYRcQBXA7NFpB75v3Yl9RjWt7VexpgE+7VfgfXPHT/72wXMN8ZcmE/cS4ArRSQMGAR8gu9rrwJAW+CVz0xgsFgnDaOB54BZdisyGfBi9YECICLXnThZBRzA+kP1dznfJ8ClItLX/qN8DKtP9Y/CAjLG7ACWAiNFJNz+R3F5Aet7gNnA/7D6pn+0y73ANOBVEWlgx99ERC4uYPefAHeJSAcRicTqaz6xn5LUl3Pbd4BX7JOQTvtEYgRWF9DlInKxXV7NPiHa1E9Vi4HDYp14rW6v31FEetrx3CoiMfb+DtrbZOLnvSylGlj/zA6KSF1gRK7libn29TXQTkRuE5Ew+6en/TqHi3VNei1jjBs4jP/PlColTeCVzztYJ+oWANuwWokPAtjdEGOB3+2vvWdgXY64SERSsU54PWyM2Za7UmPM38CtwBtYrc/LsS5rzChiXDdj9ffux0oOHxSy/gzgAuDTXF0YT2KdPP3L/qr/E1bL0S9jzLfAeOBXe7s/7UXHS1JfLo8Da7BOsu7HOnnqMMbsAq7E6gpJxmqtPoGfvzdjXft+OVaf8zas1/YtrJOMYJ2EXWe/P69jncw9ls97WRqvYZ30TQH+wupKyul14Fr7CpXxxpgjwEXAjVjfzhLIPnkMcBuw3X5N+2N9dlSAiTE6oIOqOsS69HEtEJHrH4NSIUdb4KrSE5Gr7K/1dbBaiV9p8laVgSZwVRXcj9WVsQWrL/aB4IajVGBoF4pSSoUobYErpVSIKtfrwOvXr29atmxZnrtUSqmQt2zZshRjTEzu8nJN4C1btmTp0qWFr6iUUiqLiPi9c1m7UJRSKkRpAldKqRClCVwppUKUJnCllApRmsCVUipEaQJXSqkQpQlcKaVCVKEJXETeEWtA17U5yurag5lusn/XKdswlVKqYsrIyGDVqtxjepSPorTA38N6JnFOTwE/26OI/GzPK6VUlTN37ly++OILdu7cWe77LjSBG2MWYD2sPqcrgfft6feB/wY4LqWUqpC+mP4dqxZtzJo/fPgwYLXEy1tJ+8AbGmPiAezfDfJbUUT6ichSEVmanJxcwt0ppVTFsGrzX3zx7cd4M714vV62bt0atFjK/CSmMWaqMaaHMaZHTEyeZ7EopVRImv3BXHbt2pU1b43znc3tdpORkUFcXBzLli0rkxhKmsATRaQRgP07KXAhKaVUxZeQGE/O8RRyJ/CxY8cyc+ZMAObPn18mMZQ0gc8B7rCn7wC+DEw4SilVcXgzvbz/5myS43OfBoT9x/awYfVmnzJjDAsWLGD27NkAbNu2Lau8LBT6OFkRmQmcD9QXkd1YI4q/AHwiIvcAO4HryiQ6pZQKoiUL17AtaS3TpmzH5Qjn/kF3+SzfvMk3gf/444/88ccf5RZfoQncGHNTPov6BjgWpZQqd7/P+pTqUdXpdtlleDMzmfPmVPpefxU1Gsbi9XoByCCVDC/Mnv6Nz7aH0w5kTc+YMQOPp3zHytY7MZVSVdqPG9Yxxx5oZumcOazcl8iHkz4EQPDt1961f6PPvDvzeNZ0Qcn7yJEjZZLcNYErpSqVY6kHOBC/xads3fwvif9nBQC7Vq5g1PDR7Fix3GedjKNHSU5MASDJpFuFvvm7VA4cOFD4SsWkCVwpVam88cKbvD7lw6z5TI+HT39dwbTpXwOw4IcFeB2ZLPxhgc92L78wniUJe6wZsbpOcl9ZUhqBrOuEch0TUymlytpRlztrev7MSfz6dyIAXkdmgdsddx7PU/bPxi1+1iwZTeBKKZWP9fPnkJ6W6lO2aVP+d38fSMtg5rjx+S73er1sSVgdsPjKgiZwpVSl8MmvywtdZ9TwMdQkDBywT9LYl5qW77r/G/16IMMrkxa49oErpUJKyq5NPDd0LAtmvFHwejs25CnzOjwcdKQXaT/p5lCJ4suPJnClVJW35JvPyHC5WbS+4AT7+Xuzyymi4NEErpQKSUddHg4n78l3ubds7l4vMW2BK6UqvVmvPs2vX7xZtHUnTMl3WQXL33oVilKqcpn74Uuk7ErhsnsfpG6DZgBsOBTBhlVJ/Ouq/LbKTs0erzBhaBwNaudNjoc9DtJzXFJYGWkLXCkVNIu3pLE1I5KZ46eVaPtEh4MUF6xPzdvermjJW1vgSqlKIyVhW9b0USL44PkhXNHvkTzr7Vj1Ox/M/pWz2kSTuHcffx93Zi+UitZRUr40gSulgmL6hA+yMlCay83W45F8/NpEcIYB8MOH42javJ11fbcTFmw7CDjzr7AK0gSulAqKI+IEfJ/Ql2HCsqb/2JIKWwq/Oacq0wSulCpz33/zLn8u2UHziDTufvolIKAP+quy9CSmUqrE3n3+Cd4c9kyB6yxe/D1b/rJGrjmQWqPAdfc7Na0XhyZwpVSJ7TgeRZIzvMB1fpyzlCS7X/u4ONi5+x8APE4/AxxU8ZOSxaVdKEqpMuXOcTlfhsvNx5M/IZrj4NITkqWlCVwpVSKfvvVi1vRHE0YQFhZOwq4M9ruEuLi4fLdLc3lI06tJAkK7UJRSRfLz19OZPn5o1vy63dlP9duUIqyPd7PfZfVhTxv5BBnH8w6QoAJLW+BKKR/H0tPZl7ibJi3b+pQvXLIVxEVcXByuzLACL8neY6L47IOXuOm+YWUcbdWmCVwp5eOdsaNJcoVzyTldWDZvMYKXNG8khHmz1vE4C79NPWWHh8WLvy/LUKs8TeBKKR9JLuuqkh1/rybxxBUmTm8BW/i3zyXMnftnIENTuWgfuFKKpfO+AmDt4nlZZcfSjgYpGlVU2gJXqop7K+5RdlOT339cxjFcYN/NnnbI6KNHKjhN4EpVcanuWhBmOBAGOZ9NEu+MDFpMqmi0C0WpKubnrz4m9VBgB+xVwaEtcKWqgHdeeJzk1Dq0P7k2Kzcn8/fiVQwY/Twf/G8IB8O0pR2qNIErVQXEp9bB7XKzd+MecIWT5Iyw75bU5B3KStWFIiKDRWSdiKwVkZkiUi1QgSmlAsfjyASyLxFUlUOJW+Ai0gR4CDjFGJMuIp8ANwLvBSg2pVQhPnpzLOkHUrj32Vf9Lp8xYTT/pGTq2a5KqrRdKC6guoi4sb6L7S19SEqpotqU5AZqERcXRyNvKvePGgdA6qFDzJ72IttTtcVdmZX4/7IxZg8wDtgJxAOHjDE/5F5PRPqJyFIRWZqcnFzySJVSfPrWy3z3+Yd+l8U7onn92Tj27NrO9HHPafKuAkqcwEWkDnAlcBLQGIgSkVtzr2eMmWqM6WGM6RETE1PySJVSrNt9hL9Wb+HDCSNZteyvPMsPhMG0t9/jqDc6CNGpgjgcge/HKk2NFwDbjDHJxhg38DnQOzBhKaVy27Rxfdb0lhTDl1/+mO+6R8IyyyMkVUS9e/cmOjrw/1RL0we+EzhDRCKBdKAvsDQgUSmlsmxcvZKF//cue7x1fMq9jtBP0jXdDg6HFf9BWQAOrzNkXoPOnTuXSb2l6QNfBMwGlgNr7LqmBigupaq0918fybTnHwdg1uw5eZJ35VHyQYwbSeGPtK3sStUpY4wZYYw52RjT0RhzmzFGh+BQqpTef30k2w4Y9hyPJi4uDuMoWQs1FDileMfWPiK7xV2thC334hKTN022js2/RV3QcHKBpleHKhUkU4Y/ypfvT8hTfvzAgSBEExwOipGEjXDT06PLLph8NKnbBoAbr7oDAIcJ46RWzYtVh0jJv2kURG+lVyoIPho/knhHTRK27OdK4N1XhuM+dARnVPVgh1Zhid3dEuEJ47jLTc0a4XCg7Fvhdw+6kWPHMjh6OC2rzBRx25YtWwJQr169wAeGtsCVCgpPRgYAxuFlatzD7DjsYK/UYldaOHupHeToKpZwj/2AcjtrPjFiCGe3qMmVDw8vcZ0Xn9qxyOs6nA4io4r2lJC+ffv6zN9www3ceeeduFxl01bWFrhS5ei7T99j5/rlRETVyCrbS2U9QRkYuVuZrrAwLrjrUZ+yuLg4n77nSE8YaS63lfTt3otqnjCOufyf+HRmush0enzKbrnuLhyO4nV9dOrUyWe+evWy/UalLXClyslfP37NX+u2s9fU5fjRI8EOp0IwwCXnnhzwent3juGZJ5/wewISwOn0HWro7huuACAqMyKrrO2pLWjdoXh93eVNE7hS5WTV79k33uw1dYMYSfDUcmenHGemi1O7xHB6nxuLuHXxWsPh1aN85k9tHJY13ePyyzkpolbWfLXoGlzZsyf9Bt5ZtMqN/17wsjpZmR9N4EqVkcP79pGUYD3fbdrowcRX8q6ShpkZWdO5L6Xr2akhLaqnMXhsdr/1sNFD6XvtIL91VfNkJ9uinjD0p6bdAn/ikYFcPmBIVrkjzMUdTw+md7OWiNdJjYYN6XrppdSKbUREZgSxEpVflRWK9oErVUYmvDqZDJeb87qdzJ7MWoVvEGJiM4+R4Mw+uRdV1wv5jNR26TUPFKvucJPJsVxlYvJp3ebTTQJw3yP3sWv9EqJq+38O00X33MlFucqeHv10keMMNm2BKxUg876aSWL8LgD+7503yLBPmM1fvjGYYZWZ/qNfoLGxLq2LzTzG7YNf8LveRRd2K3bdkY7sewKr4aUpmVxyet7+6P63XcOge2/KU36iKyO6XkM6nHNZsfcfKrQFrlQALJn3LfOW/c2avzbQoWcbdm5NAVf59of6E+l2kRbmKXzFE4zQ1HGQ3Sb7G0NttxAVfoRkd52sf0on9Bv5UqFV9j7rCp/5VhFp7D1a8DeSmx59iC/eeIPtXhfVxMO9cf5v4Ilt3clveVWhLXClimjHlvW8MXQYn06zktYLQ8fyv6FjANi5ZS0A+8Ic/LZyK/uDlLx7d2lDfU/2zS11wvcVuo14s9PA9Zf35d4ReUf3uW/EOJyl6YzO4fanX+KpMc/6XRbtduLwOqlVpyHO4P//81GrXg2cphpn9+yTp1++WrXgjCapCVypIvruw0nscznZuf04H7z2LMdcbo66PCyZ9y3x2w4HNbZG3sNcfvH5XPTfW6kWkX2JYpsevSG/vmNbY8fBrOlTepxdqjhquR008mb3Xjsyi/cl//Gxwxg+ahgAF958AzXcLi668vxSxRQo4eFhDBv5FH0uOzPPsvK++uQE7UJRqohMZiQ4wSMOjh85AliXAn4zbxG4nAVvHEAdW9Zh7Xbf56X0ueEO2p56GgD3DnuVkSNGYcRLRER1GpkjxBOddRVeE8ch9nitLozG3qPENG/Int3pAYkt51Umt91xDRHhJR/1PrZ1Jx4bG9guknNatqJBk4alrqdaNWu0o3pRTbnqpn8za9asUtdZEtoCVyofG1b8yfxvPs6aT3Bad9Wlu9zsP176JFBS1975cJ6yE8n7hKbVDlHL7aBr7z7cP2occSPjspbdN/xVqtuX6dVvWof/3vtksfbfPCyVtvUL709pfVInmjZpXay6y1rfO2+n04UXl7qeHud0olv7s7nr/ptp2rRpVnl5t8S1Ba5UPr787BeOudz88edYIkwmZF+aTHo+t2SXlUYcxGSG43AWbb/3PP1ygct7dGvM5hUbuOS2wpN37pwUFhnOLYNGFimOoqrjFmpXPxqw+so6kTocDq646YKs+S5durBw4UIiIqw7OZs3b87OnTvLNAbQBK6UX1OGD+aYy+pmOO5yU94Pum/oScdLGMn2Q5Duj3vNZ7kYB6aYz9LOqe/Vd9H36uz5FhGH8eRzkJLPXYeB9PDYEWW+j7LUp08fzj33XMLCrP/yt99+O5mZZT9akCZwpfyIdwT3xhtH+HEcHjdQ0+/ylrWPk77fTYKz5H3MOd319Cu+BUZADA08bnpdbJ3YFCn7RB6qRCQreQO4XK4yewJhTprAlcpl5PBRZXZ2qL4nk3DnYfaK/9vqa7idHAnLxBHmwuvJv4V9x+CxAEwd+UiRHxHi74l7+fn3WZ04tC+Bi28aULTKKxgp5lMEQ5UmcKVy+PD14RhH2Z3bHzTGuiElv2G3Tu7cnN1r13Dzg3FMH+f/Wumc+o14rdB1TrjootOJ37GpSOuecdHVecpqOY6QRnXqNWtZ5H0Gy+mX3xbsEMqFJnBV5U2OG0TLTqfx9/I9VsM7rLAtAqex6zB7PdndJJfecAfcYE3f+vhY3np+HFE10vLZunh6nXMRnJP7yR9Fd+czI1m++CfOPLdi3ppeO9PFQaeHGp4wnGHhwQ6nXGgCV1Xax1OfJ4H6JKzZYyfuwLa+G8l+8Arx+XSZ9Bv6CklJibz55iRquH2vJY+KjubhsXEBjac0IqpVq7DJG6BTmxos3FZ1xhMFTeCqijsYvweoH9A6o90uUu3nj9w/YjwAo4eNybf/uUGDhnRv35z2p50e0DhU5acJXFUpk+MGss/TmBhHCnsdtQh08o7xuLnjyUcZ94rvVR29TmtOws7s/ufGzoNkZmQPt3X5TXcHNA5VNWgCV1XGlOGPkiAx4HKzl9JfJhibmUb/0daDrU6clBw4ZqzfdS+69naf+X7Din7yURWRVL0by6veEasq56evPmTPri3EO2oWd1SufDWNOMRNjzyRY/4wjbw6zmUwde3zH5yZLs7p2bTwlSsJbYGrSm1K3CDiqc/GRf+U6IFTuUcybxl9DGO83PWE7yNX7811I0x1T1i5325f1dVt2oZho4cGO4xypQlcVUozp71IxrFjxNt93CklfFpgXWcie6lrta4dx7nz8QlF2u66W68kadf2Eu1TqaLSBK4qnZ07N/P3ntI/HrW6J4x+Y8aXaNtWJ3ek1ckdSx2DUgXRBK4qhYlDnyXZZd2B08CTAa5S3shh4Ml8Ro1RqqLQk5iqUjiRvAGSipm8YzPT6HpSbRp4MnBlhtHY7Kd5ZGqgQ1Qq4ErVAheR2sBbQEfAAHcbY/4MRGBKFdUbQ4eDq2RtkUZmH/ePfiPAESlVPkrbhfI68J0x5loRCQcC82xLpQrw3ZcfEL9tE8eOHMDhdrLP5f82daUquxIncBGpCZwL3AlgjMkAMgITllL5W7x0J16nE6hf6k5AR2TVeOiRqpxK0wJvBSQD74rIacAy4GFjjM+4SCLSD+gH1jBDSpXGu+Pj8JZw/OC7br+BtUsXsHP1Vs649EKqRUXRoXPeEcaVChWlSeAuoBvwoDFmkYi8DjwFDMu5kjFmKjAVoEePHjqkhyqVY0nHwFWtWNvUcDupVesgLVp1oEWrDnB9GQWnVDkrTQLfDew2xiyy52djJXClAurz6RPpfPr5LPh2FonFTN4NXft4IE5PUqrKqcQ9iMaYBGCXiLS3i/oC6wMSlVK25IQ9rN6czPQZn7LzQOEf1zpuGDSgH67MchyVQakgKe1VKA8CH9lXoGwF7ip9SEpli9+zrcjrNvYe4tSLL6R+g8bEOBKJpy7NTulahtEpFVylSuDGmJVAjwDFolSWn+bOYvOSBRhxAPUKXT/G46bfmOwHTN0/smS3wCsVSvRWelUhrfr9H46ExVi3hxWgkfcwHS+4gPYndymfwJSqQDSBqwpnctwgjoQVPFKOM9NFjDOB+0cV7emASlVGmsBVhbB60QL+/OYL3N4oUlwFJ+/GZj/9RmsXiVKawFWF8NOceRwOq1XodVENJYV+cdrqVgo0gasKYPWiBRwO8xa6XiPvEe0yUSoHfZysCro/vv66wOWxYQdp4Mng6gGPllNESoUGbYGroBk/dARe4+BgWMEPsez/rI7grpQ/msBV0Ox3CQVdJxgryYinevkFpFSI0QSugmL8syMgTApcp/+IieUUjVKhSfvAVbl768Uh7M8necfFxRFLMo3MvnKOSqnQoy1wVa6S9u5id7r/Pu+GJAHQP05b3koVhSZwVW4mxw0ggQZ5yh1eJw0cifSPezMIUSkVujSBq3LxyrOjOByWN3nXcDt5bOwwP1sopQqjfeCqzB3Yty/fG3WqO9LKORqlKg9tgasyM2nEQBIlxu+yGI8bh+sAF9zUv5yjUqry0ASuykyyiQXJzFNe3RPGzYMfpE69wp/zrZTKnyZwVSYmxw3E68jb+m7gyWDAmLhyj0epykgTuAq4N0c/SBL+u04GjHmunKNRqvLSk5gq4JIy/XeNODO1vaBUIOlflAqoyc8/BtTwKWscdYzwsHDqNSh4oAalVPFoAlcBM3rYGDKdvsm7S7sm/Pfm+4IUkVKVm3ahqICY9NyjZDo9PmUOr1OTt1JlSFvgqtQmD3uCRGdNn7JIj4u+l50TpIiUqho0gatSmTxiIAnOvFecDBkzNAjRKFW1aAJXJfZm3ACSJO/zTep5Ch/fUilVeprAVYk5vRF5zqLExcUFJRalqiI9ialKZPKwJ4h31PIp69zK/807SqmyoS1wVWxTXhhCgjPKp0xb3kqVP22Bq2JZt2YZ8cd8R9SJNSlBikapqk0TuCqWTz/7yme+UWYq1z8YF5xglKriSt2FIiJOYCmwxxhzWelDUhXVlGGPQY47LWMz07l/9LggRqRU1RaIFvjDwIYA1KMqsCnDHyY+R/Ku5XbQf/SLQYxIKVWqBC4iTYFLgbcCE46qiPbu2k68o45P2dmXnB2kaJRSJ5S2Bf4aMATI984NEeknIktFZGlycnIpd6eCYerb7/nMN63lpedZfYITjFIqS4kTuIhcBiQZY5YVtJ4xZqoxpocxpkdMjF4nHGomDX3aZ75d0xrcO3hUkKJRSuVUmhb4WcAVIrId+BjoIyLTAxKVqhBmvPUyia6IrPnq7jBuvvexIEaklMqpxAncGPO0MaapMaYlcCPwizHm1oBFpoIufluaz/yTY58NUiRKKX/0OnDl15txAzkSlj2ivN5pqVTFE5Bb6Y0x84B5gahLBd/Pcz/1GZQ4RvTks1IVkT4LReWxcPG6rGlXZhgDR08MYjRKqfxoF4ryMXnEIJ/5u/vdEqRIlFKF0QSuskwZ/jAJkj1yfKxJpnGzlsELSClVIO1CUQC8/Uqcz92WDT3p9B+jXSdKVWTaAlcAJO8P85l/YIw+50Spik4TuGLet3M45nJnzV97lT5UUqlQoAlcMW/R8qzpDi3q0fG0HkGMRilVVJrAq7hJQ5/ymb/hrgeDFIlSqrg0gVdhE0YOItFVLWt+0IABQYxGKVVcmsCrsBSTfclgtw4tqZ0jAsIAABuqSURBVN+gQRCjUUoVlybwKuqloWN85q+44c7gBKKUKjFN4FXQPxvXkebyZM3ffON1QYxGKVVSmsCroBkff5o13alVA9qdfGoQo1FKlZQm8Cpm4sjsZ53UdRuuuV1PXCoVqjSBVyEpSUkk5zhx+dDYkUGMRilVWprAq5CJEydnTcdmphWwplIqFGgCryLmff0pRrwAxHgy6D/6pSBHpEJdUsJeUpISgx1GlaYJvApISUxk3tLsQRoGjnkuiNGoUPDIK7PpN+GrPOUbV63k+vd+4fcff6b72kTOWb7b7/Zfz/qOY0fTfcpeeXkmsb+uZMbb/1ekGNxuNykph3zKnnrjS2J/Xcnn/7ewwG137khk3ZrtPmXvz5zP/PlrCtxu5fo5pKUfynf57CU7ue6bVUz5bXOB9aRlevk2+SBurylwvdLSx8lWAR+NnwT2wwZjZH9wg1FFdujgAY5nZNCgQcNC1x3ywkeEe48z5pm7i1R3eprVhVY9MpKl8xeyd28iG3cc5JVePei59zBLurYB4H8H9vPDnJ/ZkHCMptGwPj2CBd3bcSQxFXdD4YDTyZFDh3C5XKz4czUfrU/ms07NoUEs17z9IxMfuoIPp37Bye2bMbd5CwDmeiOJ+W4JP67fwwddW9I0NYNrt+ylSd0aPNGsHgARmYbjTgHg8n8SqJdxnKdv78t7Ha06BtSqwddTv2FDoxpEHXPx+ikNSfCsZv3az5jvvY/fmtYC4IdVW+h8WmsAnoytBd5MRsz4hXUe4fOmtdl8RhsMHpJ2HGLrnk3gGMyXKxuwteG7dKwTzZcJ+9kSvpGn0+azMuohXnOkQyQsdKdyxf40qtUIJ1yEKJeTcStnUqN6I5xRpzJuWwIHPdljyv7e62RaR2bf9RwoYkzZ/ofIqUePHmbp0qXltj8Fa1csYvaX32bN6+DE5evg/n0MmfEHZ8le7hh4f5G2eXbsO6SFR/F7+5PYWSOchH914bmxb/JT+y7UPXqc5Y1rEZFpGPD3Era66xKVkcpbZ3YDoP8fSzgSXoOfTmlDYqSL17avYnUCtIg8xojOvQrcb0y6h+Tqvm26q1dv5/POLUt07BXFqWY166RzvssvMN/xk/wbADFeIjjOMake8DgS/tWlxNuKyDJjTJ6nzGkCr+RyJmxN3uVr5aJFTPtzL5+ddhIAfbamAIZXLzmVxQt+474G7ThlXxoxR9JpHb+LT3p0plFqBpvqBL6lpoJvTe92xERElmjb/BK4dqFUYhPiBoI9unysSQluMJXMtAlTWO2uz6ddWtMlMZXWCcm4w1zMOaUZAFeu3cmXHZuDnbwBfmllXcJ52sYkaNAOgPX1IqFeJPNbWl0HmrwrrxX7dnFR4/YBrVMTeCX1xQdvkmIn77oeQ/8xE4IcUeg4fOgQ93yxjI47N7KmRQd+a24NNXfujv0saFHXWunU7O6IlQ2jWdkw2qeOLzs2L7d4VWjoWiMi4HVqAq+kVm1Nypq+9r47ghhJxbZ182Zem7OKk81eNnkbUD0jDbcrgoWnn8LCFr191s1K3ipf9Uwy+yQm2GFUSDXCwgNepybwSmjK8IfAYSWbRt5DNG7RKsgRBdfIkRM5XL0uLw+5ic8/mkVSYgoft+jGzpoRpLkc0LU10DrYYRZbhDnGccnb5dLa/MMWaZc17zRu7mIayzidFZLdjTrEjGYXLUglmq/k6jz1XGi+ZRG9SacawxjOfuqyiLP4U84G4BHzEqeyhkjS8OBkLlfwb77hLmZm1XGtmUFd9rOI3qySbrQ2/9COjXwrV/iNFeBm8x7RpDJVBnGPmcwW2nAFn/OovJnva9HU7OQQtbiDt0kjknekP+3MBv6RDgDcYabhJiyrXoDp5hqSachKuvEV/2W/WF1cw82znMRW7pKZPvvoZX7nYubyB+fgwckNzGAvTRgtY4gzT1GfZGpzkJ+5mHelHxebb/heLiXGJPIaA4CCL30sCT2JWQnpiUtb+gEW/Lac6131ym2XUeYIR6VG1nxns5wnGct8/kVzdtCC7XzCzXRlCaPkOZ4xI2jPBnbSklZsYQut2UtTzmY+BuEAdZnEQ7RhE//hK2pxCAO4CSMcN14cLOB8ltKLm/iAJuzBAAeoS118Lxk1wG6a4SaceiRTi8MkJ0QTE5vKLfIZF5tvuIZZ/MC/2cFJPMI4VmzoyCmtbiciYkhWPSs2nE7XDosByEh3cujYYGLqjMtanpIUxQY5F++uUzmvm1W+ZtO5dGq7IGudo4fe41DqcS67rCvTPnmUVk3WkBLZkDrsZ//2ThjvNnr3Hs2fi0fgldq0arCKZxxTueqfFFrF1GTp8TUkyCm0OL6W+s51dG75By6sy/bWxl/A840f4CHzP3rxF6kHI4iufRyA3XvrExFpSK8dRXN20rLNRyz+60WiwzayPW0iF3duwa7kiwBYQi861HqB3o0OEl3vFA4mLGX9llsAcB+MJab+FbQ55TqOHP+L9PSdZGSk4PUe59RTXiYtbRu790wnJuZiVq26G2MyOav3b1Sr1qhEnyu9CqWKmDT0SRJd1iVQzaI93PP4mEK2qDx6zllCr627mZ3jxGFZed/cwHJ60IAE6rCfWhz2WZ5BGJN4iBv4iFgSilzv0YPh7Nlfj3at4rPKlq7uTlRYYzp0sG6sWbbmdIY8PJNNG1awePnd7E45mQfvfYvIqChee2MgLRqdylXXDmDdqiX8uugJGtdJJsPdnz0HvsPpzOCR/j/y+adTufq6fn5jeH3ykzgci7n3trlUr259lt6e/gJHUn/h4j6v0KFdR6Z9MJpDqWt4fMAnAKTsS+Kjz4bR/47xRET47+t9690zaNpwH42bz6Rzx8COu/rzL9Y3qDotPqRb697M/PVlzu1yPU3qNMta54cVs+nV/kJqRdZixsy+NGy4nb59tuSp6++t8zkUv4L27a+gVn3fb6/LF7xEVI2mtO96c5Fj+/ufUeze/T7nnbsSl6tG4Rv4oQm8Cpg+6QU2Jx4DoJH3APePej3IEZWPkSPfYNK555S6nh5mEdVJ4zRWUJ00urAia9kEBhNLPFfxCYLBQdH+brbvrEfL5vsAOHLwCTbt+Ya61buTcnQVXq+LsLA06oR15yi/8Z/zXqXtyadlbTtx2rk0qZvCf69ZD8DcL98m03i5/L/3lfpYK5uNe9bgdLpoG9sh2KHk4fV68HgOEh5ev/CV86EJvAqoCl0nM9+axtgm3eixM4nv2hX/6+iNy9eA00uDg/HEd/Syo25tbmI6HlzEEo8Db5HrWrzyDCKj4snIiEQcmUSH9eRA6nZczpqc1Kwz117tv4WrVHHpdeCV3KQRA0CsMS0bRx4LcjSB9dRzH1LjeDJ76rbi8849AYqcvB9Y8DsjRgzkhf/dQ8/u86Br8fe/amMbGtZKpHWrcURGRtP5tDMA6Nun+HUpFUiawCuJRDt5h3nC6DckLrjBBMjNb/9k3fxyZqdC12198DhbakcQt2wB/R9/KKv8hdUj+fmXV+jZvfD9ZWY4cIZnt8C7dP6DevUbaqJWFVaJE7iINAM+AGIBLzDVGFM1Ol0rmLgRcWA994f7H7w3qLGU1pRXX2ONtLBORLYqvM/w6tXbef7286hVpw4b167kp5qf8eL4X+jR0XrqXM8Ccv+ydafQuM5ukg43oH3jm7n0v3q9vAotpWmBe4DHjDHLRaQGsExEfjTGrA9QbKoI3ozL7jppQBL1Gxb+5LqKZsnC+Xz9y3qmnHsmdDm/wHV77j3Mkw3TOPsC6+FD/KsLH30wnlT3TNqclESHNgD+H3F6wrpNzXno/l+1Za1CXokTuDEmHoi3p4+IyAagCaAJvBwl0SBrekBc/jc6VFTdv1rCnuhacO6ZfpffsmQjD13Xk79XrOCiq6/NKk9K2MtbHz5Dr+4LiW1a8D6MF9q3/5GaNerizkinb5+SXYurVEUTkD5wEWmJdXpokZ9l/YB+AM2b6/MhAmnKsMfAaV1XGkonLj+cNJH363dlbf1IiA7zu84NyzdR70g8w+Os/uwWJ1nX+W7bvpn1qy+lWrSHXoX0ay9afhYD7nmd2nXq5CitHYhDUKpCKHUCF5Fo4DPgEWPM4dzLjTFTgalgXUZY2v0pS3JCAvF28nZmuug35IUgR1S4uyfNZe7JjeHks/Jd59f6hg6dukKuZyd/O3cW4dWeAaBatL8tYcX69pzS8nYuuexGQK8SUZVfqRK4iIRhJe+PjDGfByYkVRRTJ7yd9e71H1ixT1y+M/51nul0Hpzc2O/yd5P/5j/X35CnPGHvblav/hfOcC/hBTxltXmzr2nbtoMmbFXllOYqFAHeBjYYY14JXEiqMF/OmITb5QaggewjJjY2yBH5N/uDdxjUrBt0Oi/PsivX7qRlykaeHjkEyG5tH9i/n09mvgHVvqfNSYk4/TzA7WBKdfpcsIAPP3yJ/1xyB23bVry775QqDyW+E1NEzsZ6vNYayLp97RljzNz8ttE7MUtvz47NTHt3OgAxHg8Dx1SsZ50s/Gku3/0Zz9tn+++gfvL37xg89Cm/y158cSA9en6Xb91JCdH86/xviG1cyFlLpSqZgN+JaYz5jayrj1V5efvtmeCwpq8fNCi4weTy5rj/Mar7hXB23q6Sh379iWdGPZ6nbxvg3bdfpvlJb9Kjp/9616y4kkce0y95SuWmd2KGkHdfj8NrJ+9Yk1Jhuk7GjhjHG+dfAN0vzLNsSvJarrz+1jyJe8vWTcz6bAS9ui+iuZ+HBx5IjqRx7Cv0Pu9C7dtWKh+awEPIjgPZ0/1HVowh0mJ/XQnnX5Cn/Inf5vLYsGfI2b8NsH/fPlasOh3A72WAS1Z15anBs8siVKUqHU3gIWLi0KHgst6uru2Cf7dl7K8r85SNXPkT9w9+3JrJ1eIeN24wXbvNybe+1INxXHn1bdraVqoYNIGHiGQ7eTfIPM6VNz8QtDgumbmQ5bF5H0r/9oGNXHoiedv2p6Qw+e3H6NXzN7p2y1vXzl11ueuOJWUVqlKVnibwEPDqs6PAvmFxwOjngxLD6BGvMPH8PpAreT+75HseHPIkubtKfvqpNeKAXn5OTK7f3IwH+80rs1iVqio0gVdwX86YxKEw6yrNGJLLff9xcROZfN5ZcL5v38ZD877nmZFP+nSVvPTSIE5p9xPVa7oRR966dm4bwF33PKbdJEoFiCbwCm7FP4kAiNfJwFETy22/n8+azoAGHeE839ve+y34k1EjHvBJ3K9OupDO7bfS3c8wh4uWnEOrRt258fYHyzpkpaocTeAV2JThj4HD6rK46qq+5bLPfcnJjPpgHrO6dfQpv+uPldx93em0HWH1v7/w6jVER6XQoc1uOrfPW8+W9bdzxX/vom8ffYCZUmVFE3gFtWzRL8Tbybs+SXTu2rtM95eSlMSEN2cx+bxzoFtbn2ULm4TT9tk7ARg77ibO6LaYnqflrSMj3cmhfQ9y4+0PajeJUuVAE3gF9fOcPyAM6rm9DBpbts/5/u2X77lWGsJ5viO7f5C+g4suuZJ9yclMfuss2rZK4Aw/V5MsWdOZpx7+okxjVErlpQm8Apo0YgBpYdZADTc+OKBM9/X4izOZfrrvw6AG/PoL//l3V7YnDeXnXx4FoG0r3+1SD0Rw9rnzqRcTo61tpYJEE3gFM+/b2VkDFMdmHi2z2+Xfn/wmT7bvDbmS91N/P0qn83dw+BjUjcm73V9LzqH/3f+jXoyfhUqpcqUJvIJZ+MdGcFrT/Uf/r0z2ceu0H/ipvW+f+jRzK5GkQ7u86y9dchFPPjkJ0EESlKpINIFXIBPiBpHptEZi735K4MdtfGvCaww99Xxokz2O5stmELHW0KY+du6qyxWXzNUuEqUqME3gFcTGtctIwUrejbwHufz6uIDWf+aXf7Ht1POz5luaLYxliM86K9a3p22Da7jy+nsCum+lVNnQBF5BfDz7q6zp+0e9FrB6p717BsNaToaa2WOSTTfXZD3IfcWGdjw+8FtAu0eUCjWawCuA14eOzHonzunetuCVi+C5F2+ne8+/eJmnWd1yclb57eYtLuZblqzpRHRmOwY9+pImbaVCmCbwIHvn5aEcsJ80WF9S6Ht5XInqSU5MZNr0h+jVdSkn96zFHfKJz/KBi1/nqksG0LHTBE3aSlUSmsCDaNWy39l5xHoLaridDBpb/EEaFi6YS1LCEGrXT+f0rnCLfOazfMC8eQwf+Qj8692AxKyUqjj8PDNOlZcvvvoxa/qxscOKtW1SQgJvv3c6GZ4HqV0/nb005tYcyfvsnQeY38BhJW+lVKWkLfAgmTB0GLisC74vvfDsYm370vgr6N5xHS2bgxdhDlfzqdyctfz++b8xMq5iDXislAo8TeBB8OqzozgUZiXvhq6D9Dwr75iS/jw37kZ6dVtCd/tBgd9wOTPkzqzl163cyhuDr/Y78rtSqvLRBF7OJo8YxKEw63rvhiaZB4YW/ozvqZPH0Lrdu/SyHySVSEMeFd8HXL2VtIHLBt8U8HiVUhWXJvByNGXEgySIlbzrmyQeGFnwUwaTEhJYs/4sWtu3t3tx8FLqa6yp0SRrnRtWbOb1R68l95BmSqnKTxN4ORk1fDReRz0AYkwiA0dOKnD9lydeTJcOm7PmfzEX8rajP+QYkjLhX120u0SpKkwTeBmbNGIgiRKTdb1P48h0+g3JP3k/9/L19Oq6jC72QwK30pph8hJZt04Cb+/bzKXXXluGUSulQoEm8DKyYukCfv5iAalh2Y9djY1Mo9+Ql/yuP/alWzmjx5/06gppVOcFhrNFfB8N+Oj8rxgSNwztLlFKgSbwgMtqcQOEWb8aeDIYMOY5v+s/9/IN9Oq6lDN6wCbaESfP51nnqT/m8sizz2h3iVLKhybwUho/ahDhmYaEE0lbslvcrkwX553XiXMuuNJnmxNJewtteL/rU4yXOn7rHrb0ewY+8aQmbqWUX5rACxC/ezdzZk3CHDlApqcOya4wP2vV9+mfBqjjhofHxgGQsHcvo0bfTZuO61lc6yyW04N/uj2d7z77LVzEPbf2oUWrtpq4lVIFEmNMyTcW+TfwOtYYMm8ZY14oaP0ePXqYpUuXlnh/ReH1GhL27OarmRNxpx0mzFsNY8JIcFbHmeki0ms4EpZZrDoN4HE4SaxZB68I1TwZpETV5kh0dRwRGUR7U5nfoOSjxt+8dCOvPHFjibdXSlVuIrLMGNMjd3mJW+Ai4gQmAhcCu4ElIjLHGLO+5GEWzTdffMKqP3/DEVkdb3gExyMiyIhwcjjawfHqXtwuF/SowfdyI7U9Bznoql1gfTEmkSiTSoaEs1ea+SxzGg+ZUrovKmHmOG6JINKTya1/LKaOax+Dnx2avYK2tJVSJVCazHQ6sNkYsxVARD4GrgQCnsCvmzOVjVEtOSbVOEY13LXbwX/8DN7oR2HJG6AWB4mWVBxkshffBN6X79lv6rFUzvC7bTWTxv1Hp7Ch2im03JtEeMphure9kP9ccb3/nV3YvUhxK6VUYUqTwJsAu3LM7wZ65V5JRPoB/QCaN29eoh3VTkvjJOceIrwZRGRmUC3TTYTHTYTbQ7jbjeu4B0eGB1dmKngy4Gg4SBo1ImK5/pYH+ParD+l02pk0bd6GJs2a+dlDQS3gorSOS959opRSJVXiPnARuQ642Bhzrz1/G3C6MebB/LYpjz5wpZSqbPLrAy/N88B3g09/Q1NgbynqU0opVQylSeBLgLYicpKIhAM3AnMCE5ZSSqnClLgP3BjjEZFBwPdYlxG+Y4xZF7DIlFJKFahU18cZY+YCcwMUi1JKqWLQMTGVUipEaQJXSqkQpQlcKaVClCZwpZQKUaV6mFWxdyaSDOwo4eb1gZQAhhNMeiwVT2U5DtBjqYhKexwtjDExuQvLNYGXhogs9XcnUijSY6l4KstxgB5LRVRWx6FdKEopFaI0gSulVIgKpQQ+NdgBBJAeS8VTWY4D9FgqojI5jpDpA1dKKeUrlFrgSimlctAErpRSISokEriI/FtE/haRzSLyVLDjyU1E3hGRJBFZm6Osroj8KCKb7N917HIRkfH2sawWkW45trnDXn+TiNwRpGNpJiK/isgGEVknIg+H4vGISDURWSwiq+zjGGmXnyQii+yYZtmPQkZEIuz5zfbyljnqetou/1tELi7P48hJRJwiskJEvrbnQ/JYRGS7iKwRkZUistQuC6nPV44YaovIbBHZaP/NnFmux2KMqdA/WI+q3QK0AsKBVcApwY4rV4znAt2AtTnKXgKesqefAl60py8BvgUEOANYZJfXBbbav+vY03WCcCyNgG72dA3gH+CUUDseO55oezoMWGTH9wlwo10+GXjAnh4ATLanbwRm2dOn2J+5COAk+7PoDNLn7FFgBvC1PR+SxwJsB+rnKgupz1eOuN8H7rWnw4Ha5Xks5f4hLMELdCbwfY75p4Gngx2Xnzhb4pvA/wYa2dONgL/t6SnATbnXA24CpuQo91kviMf1JXBhKB8PEAksxxqzNQVw5f5sYT3X/kx72mWvJ7k/bznXK+djaAr8DPQBvrZjC9Vj2U7eBB5yny+gJrAN+2KQYBxLKHSh+Bs8uUmQYimOhsaYeAD7dwO7PL/jqXDHaX/17orVeg2547G7HFYCScCPWC3Og8YYj5+YsuK1lx8C6lEBjsP2GjAE8Nrz9QjdYzHADyKyTKxBzyEEP19YvQLJwLt219ZbIhJFOR5LKCRw8VMWytc+5nc8Feo4RSQa+Ax4xBhzuKBV/ZRViOMxxmQaY7pgtV5PBzoUEFOFPQ4RuQxIMsYsy1nsZ9UKfyy2s4wx3YD/AANF5NwC1q3Ix+LC6jqdZIzpChzF6jLJT8CPJRQSeKgOnpwoIo0A7N9Jdnl+x1NhjlNEwrCS90fGmM/t4pA9HmPMQWAeVr9jbRE5MRJVzpiy4rWX1wL2UzGO4yzgChHZDnyM1Y3yGqF5LBhj9tq/k4AvsP65huLnazew2xizyJ6fjZXQy+1YQiGBh+rgyXOAE2eT78DqSz5Rfrt9RvoM4JD9Net74CIRqWOftb7ILitXIiLA28AGY8wrORaF1PGISIyI1LanqwMXABuAX4Fr8zmOE8d3LfCLsTok5wA32ld2nAS0BRaXz1FYjDFPG2OaGmNaYn3+fzHG3EIIHouIRIlIjRPTWJ+LtYTY5wvAGJMA7BKR9nZRX2A95Xks5X0Co4QnCy7BuhpiC/BssOPxE99MIB5wY/03vQerz/FnYJP9u669rgAT7WNZA/TIUc/dwGb7564gHcvZWF/fVgMr7Z9LQu14gM7ACvs41gLD7fJWWElrM/ApEGGXV7PnN9vLW+Wo61n7+P4G/hPkz9r5ZF+FEnLHYse8yv5Zd+LvOdQ+Xzli6AIstT9n/4d1FUm5HYveSq+UUiEqFLpQlFJK+aEJXCmlQpQmcKWUClGawJVSKkRpAldKqRClCVwppUKUJnCllApR/w+oADJOmPZ/tAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "continue_train = True  # set to True to continue to train\n",
    "num_big_steps = 30     # number of small steps\n",
    "num_small_steps = 200  # number of big steps\n",
    "if continue_train:\n",
    "    for k in range(num_big_steps):\n",
    "        for j in range(num_small_steps):\n",
    "            dXY_list = np.append(dXY_list, np.zeros((rep, 1)), axis=1)\n",
    "            dX_list = np.append(dX_list, np.zeros((rep, 1)), axis=1)\n",
    "            dY_list = np.append(dY_list, np.zeros((rep, 1)), axis=1)\n",
    "            mi_list = np.append(mi_list, np.zeros((rep, 1)), axis=1)\n",
    "            for i in range(rep):\n",
    "                minee_mine_list[i].step_minee()\n",
    "                dXY_list[i, -1], dX_list[i, -1], dY_list[i, -1] = minee_mine_list[i].forward_minee()\n",
    "                mi_list[i, -1] = (dXY_list[i, -1]-dX_list[i, -1]-dY_list[i, -1]).copy()\n",
    "        # To show intermediate works\n",
    "        for i in range(rep):\n",
    "            plt.plot(dXY_list[i, :],label='dXY')\n",
    "            plt.plot(dX_list[i, :],label='dX')\n",
    "            plt.plot(dY_list[i, :],label='dY')\n",
    "            plt.title('Plots of divergence estimates')\n",
    "        display.clear_output(wait=True)\n",
    "        display.display(plt.gcf())\n",
    "    display.clear_output()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Continuously train the model with MINE. The following can be executed repeatedly and after loading previous results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAEICAYAAACtXxSQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd5gUVdbA4d+ZRM4MOWPENMIIIoqAmNPqmjCjLqsiIqgoGSQqBsSAi4nPnNOirgpixMAgqEiQnASGKDnMzP3+qJqhc6yOc97n4aG6wq0z1d2nb917q0qMMSillEptGYkOQCmlVPQ0mSulVBrQZK6UUmlAk7lSSqUBTeZKKZUGNJkrpVQa0GSeRkTkKxG5JcExHCkic0Vkp4jcGcL6I0TkFXu6mYjsEpHM2EeavkTkGhH5PNFxqPjSZJ5iRGSliOy1k95GEXlRRKqGWUYLETEikhWDEAcAXxljqhljJoWzoTFmtTGmqjGmOAZxpSVf76Ux5lVjzFkx2l/CKwzKN03mqelCY0xVoC1wEjAkwfG4ag78keggXMXoR0uppKLJPIUZY9YBnwLHei4TkQwRGSIiq0SkUEReEpEa9uJv7P+32zX8jiJymIh8LSJ/i8hmEXnT335F5CIR+UNEtts1taPt+V8CXYEn7XKP8LFtS3s/O0XkC6Cuy7KyWqaIXCUiBR7b9hORj+zpCiLysIists9QnhGRSvayLiKyVkTuE5ENwIv2/AEisl5E/hKRW+x9HRZGeXfbx3K9iPR0iauSiDxiH+u/ReQ7l21PFpFZ9rH6VUS6BDiujUTkXRHZJCIrXJupRKS9iBSIyA47vkcDvJc3ish3LtsaEbldRJbYx32UiLQWkR/s8t4SkRx73VoiMs2OYZs93cReNgY4zeX9fdKef5SIfCEiW0VksYhc4bLv80Rkgb3fdSJyj7+/X0XJGKP/UugfsBLobk83xaoFj7JffwXcYk/fBCwFWgFVgfeAl+1lLQADZLmU+zowGOsHviJwqp/9HwHsBs4EsrGaVZYCOZ4x+Nn+B+BRoALQGdgJvOIZF1DZXna4y7azgavs6YnAR0BtoBrwX2CcvawLUAQ8aO+nEnAOsAE4xi77ZXtfh4VR3gP233wesAeoZS9/yv67GwOZwCn2fhsDW+z1M+xjtgXI9XFcMoA5wDAgx37flgNnuxy36+zpqsDJAd7LG4HvXF4b+2+rbv/9+4EZ9j5qAAuAG+x16wD/tI9RNeBt4AOXstzeX6AKsAboab9vbYHNwDH28vXAafZ0LaBtor9D6fov4QHovzDfMCuZ7wK2A6uAp4FK9rKyL5r9Zb3dZbsjgYP2F85XAngJmAI0CbL/ocBbLq8zgHVAF88YfGzbzE6KVVzmvYaPZG6/fgUYZk8fjpXcKwOC9YPS2qWcjsAKe7oLcACo6LL8BezkbL8+zN7XYSGWt9fjeBUCJ9t//17gBB9/733YP6Au8z4rTZwe8zsAqz3mDQRetKe/AUYCdT3W8fVe3oh3Mu/k8noOcJ/L60eAiX7eszxgm8trt/cXuBL41mOb/wDD7enVwL+B6on+7qT7P21mSU3/MMbUNMY0N8bcbozZ62OdRljJvtQqrERe30+ZA7CS2s92E8pNftZzK9cYU4JVM2scQtyNsBLDbo+4/HkN6GFPX41VQ9wD5GIl9Tl288V24H/2/FKbjDH7PPa9xuW163Qo5W0xxhS5vN6DVUOui3Ums8xH/M2By0vLtMs9FWjoZ91GHusO4tD7dTPWWdEiEZktIhf4KCOQjS7Te328rgogIpVF5D92k9EOrB+RmuJ/hFFzoINH3NcADezl/8Q6M1llN691DDNuFSLtGEpff2F90UqV1oo34iPxGmM2AP8CEJFTgeki8o0xZqmPco8rfSEigtXcsy6EmNYDtUSkiktCb4ZVc/Tlc6CuiORhJfV+9vzNWAnoGGP1G/jiWeZ6oInL66Yu06GU589mYB/QGvjVY9karJr5v0IoZw3WmcDhvhYaY5YAPUQkA7gUeEdE6uD/2EXqbqyzuA7GmA32sZ+L9UOPj/2tAb42xpzpJ+7ZwMUikg3cAbyF+7FXDtGaefp6HegnVodjVWAs8KZdu9wElGC1mQIgIpeXdnQB27C+tL6GCL4FnC8iZ9hf0Lux2mBnBQvIGLMKKABGikiO/aNxYYD1i4B3gAlYbdlf2PNLgGeBx0Sknh1/YxE5O8Du3wJ6isjRIlIZq226dD+RlOe67QvAo3YHZqbdCVkBq5noQhE5255f0e5MbeKjqJ+BHWJ12lay1z9WRE6y47lWRHLt/W23tynGx3sZpWpYP2zbRaQ2MNxj+UaPfU0DjhCR60Qk2/53kn2cc8Qa817DGHMQ2IHvz5RygCbz9PUCViffN8AKrNpjHwC7qWIM8L19anwy1hDHn0RkF1ZnWV9jzArPQo0xi4FrgSewaqUXYg2VPBBiXFdjtQ9vxUoULwVZ/zWgO/C2RzPHfVgdrz/azQHTsWqUPhljPgUmATPt7X6wF+2PpDwP9wC/Y3XQbsXqeM0wxqwBLsZqLtmEVYu9Fx/fO2ONrb8Qq416BdaxfQ6rgxKsDtw/7PfncayO4H1+3stoTMTqMN4M/IjV3OTqceAye6TLJGPMTuAs4Cqss7YNHOp4BrgOWGkf01uxPjsqBsQYfTiFKn/EGk45H6jg8SOhVErSmrkqN0TkEvvUvxZW7fG/mshVutBkrsqTf2M1dyzDaru9LbHhKOUcbWZRSqk0oDVzpZRKAwkZZ163bl3TokWLROxaKaVS1pw5czYbY3J9LUtIMm/RogUFBQXBV1RKKVVGRPxeMa3NLEoplQY0mSulVBrQZK6UUmlAk7lSSqUBTeZKKZUGNJkrpVQa0GSulFJpQJO5Uipib70yhXufG5voMBSazJVSUXiubmVebn0eDz0yNNGhlHuazJVSEduWVR2AgxJkRRVzmsyVUipGpqwp5OQfF8RlX5rMlVIR81chv+OVh2kwcx4ff/R2XONJNsOW/sXKvaE+UTE6msyVUlHzfCrCtw2PBWDRikXxD6ac0mSuVBDtp3/EiTM+SXQYKWWb1HS8zDFPj+Le58f5XNb35QkMeHG84/v89JtP2bNnN4++9TR93n7c8fKdlJBb4CqVSlZnNkt0CAGtXrGEChWrUr9hw0SHAsCYicM4cMKljpf7xNEXAjDBx7I3m5wJwENR7mPqB6+zaP8Wxl95B1/+MIOexQ05a/rLfJ57CgC3LV/M3LWLuabzRVHuyXlaM1cqxbVfuZvLf/8uQXv3fuzk1ppV4rb3TRv/4q03nve7fOkf87j7pQns27s3pPLur3E0U+udCsCqbRsB+LVKi7Ll/1i5nruLD/243zlzKnfMnArAHf/7D09++ya7D4a2L6c5ksxFZKWI/C4i80REnzqhVJz9mX14QvcvLj2h8Xqu8JPPPcZxCwq5s367snkXf/wC/V86VD8fvPQnXm16JhNenRjxfjZmNCib3uHRfPQWebxDHr+vW8g7FTowuuhI2n33s1cZDWbOY+P+gxHHEAona+ZdjTF5xph8B8tUKirffzWD3q8+wltTn0l0KBF5fvLDnPLFe4wdf39I649+ZDAjJzp/Ac89L4zjvmfHeM0XH3nb1wiXo778ko5fvA9Yia3HB5Ojjml+Ve95P1Vuy2tNzyp7vTO7IgB7M/0PhL//5cdoMHMeJ8z4X8SxDFg4t2x6OzV8rvPj37siLj8U2syi0to7Kwt4t9EZfFo1trWiWFmcfYDlWa34/vijQ1r/ybaXM/mEfzoexystz+X/Djvf7/KFzRsy/tFhbN+2DfFRM98utVmR1bLs9cwaHb3W2bDxL9546zmv+Vd/MJkGM+d5zQ9W///0g5f5pdIJ9rruyXzRvF8p+P5bAKY26Qq418ABNhzwn3xH/+//ePDr1/wubzBzHhN+n+Y2b9me/fy6c0+QqCPnVDI3wOciMkdEevlaQUR6iUiBiBRs2rTJod0qFdi+CtnW/1nR9/Xf/PYk+rw8gQYz5/Hko8OjLi8cnonr1alP88vsWRGXN+DZMQywa9qPPDKc4U9E9/dMr3UKE0+8lAc++I/Xsis+OjTv6Bkz/JZxz/cfcFduPm+/O9Vt/pc+Ej+ACXDV6dDnx9OzxnFe88/59GU6ffEuXbYZLjhQzX8BwP6SYr/LnqxwAo+VtDkUi4/zkUc2N3F7/dCKDZxd8CcztuwIuN9IOTWapZMx5i8RqQd8ISKLjDHfuK5gjJkCTAHIz8+PT6NaHH32yav8UfAL/YZOQERPeJKNE1ebf1y3c9n09iL/F4IsWfQblSrVoEnz5mXzpk5+mB9qZTLgxO60PtI7yfjl45sy+uHBPNnuco7fMp/POcXnZoUbNlCvQQOfyx6aMIiX8q+wpoEnTzybvVKZkfbyCQ8PIScjk779R/rc/rFHh/NxmyMZlFHF68BurFmdhtv+dpv3TbUOZdPbMup4lTen4HsmrP2dH2ueCMCLlYTLfe7Zk/939RcfI3s6fvE+KyqGfuyfqeX7R8SXBRmtQl539b7YXETkSNYxxvxl/18IvA+0d6LcVPLS/h08dNp1PPGQ3nAo0b6fOZ0nJgwDgp+K+zN+7L00mDmPYRHUWE9bX8I98z51m/dx42p8WL8rU775yOc2r77wBDM++9BvmUUZmRw343MGTR7Fk+2sVPdbhWPd1un96iNl0xOmvQDAow8P5YzPXueVqU8A0OvNiTxqJ/JSe6Vy2fSSPxfySLvLGHfiJYx8fJhXHK2//JYHT7yE+RXa8Pmq+Xge4f1Z2WyqcajGu6Vi8J/R55b8xFc1TmafVALgl0oncO6nLzHw2cB3YyzO8J++5tjNK6W+b9zKrZknmMv/+2zI6wIckAohrxurql7U5YpIFRGpVjoNnAXMj7bcVDOrhvXh2bU7tp0cytvY0fcywOWLP37PGsbk2+Oc7WEWvjrqACaOG8yoCYNoN30aQ546VBPd3KA2AHNbtfC5nRHfSaqkxGqb/6r6yW7z92XmALCnYo7XNs8/NY67W55G36xKh+Y98RANZs5jQeNGVjxZtdiUUY8XjrrQbdurXToS3210Rtn0y63PY9GCX5l1RDP+yDmaryqX8MKzD/NRvS5u2/d6032UR5e1u8umJx9/KRvWrWaoy3HZLe5NE3szKrm9/q7qSXxW+9Sy1y8ccYHX3+vJ17GcW/F4XjzsPMY9PcrnNm++85LbmVIw4Y72+bbqSWGtH44PC7fHpFwnmlnqA++L9YZkAa8ZYyLvFk5R0bQbjR3el1/aHUvXP9fR+54RToXk03uvv8jixb8ycETkQ7WSzbOnXMJeqczp77/K2tVLmH38P8qWGR9Trl466RT+ymwMwPtHVWS0Pb/ETjD+krbnGX7+9P9Sq2gHvXcegLonui0bP+5eZp98DQBvNzmTJ4AVf87nrLVbuHD1D2ytbw3L2JxRr2ybdQesdtWCSnn+/3D8tycDDF/yI0uqWk0903JPZ1qu9zquyf2a9ydTXNO9vOHfvsuHbS72Wf7UAB2ioRg2eRRTjroQ6nf1u87jR1/oNe+EGf/j1APO9Lv56liNtQW7YjMOPepkboxZDpwQdMVyQiJI67NPOp4fqrSjWd2tMYjIMuvr//HojjWsqZvLqgY3cuTrL3Jpj55hlzPi8WFU2LaLgSMejUGUliULrLvMHd6mjdeypQsX8tgvn1Br+y56nnElD837lL12MhhbtRLLXBL52HEDoKXVdrozp5JXWUBZIgcodunrKE3i/pLpU+2uoMKYATzT8SJKyGRfZlPWZsJLBw9dZtFg5jy6/v0DP3S4xG3bn76dwca/VrGzXlvebtaNM7b+VLbspOkfsSazGT3quNeAXeMM1dfVOwRfycWMmt4/DB8GSLTRmnKUd6IOxcaMBrzb0Hd/QCrI8ldBiLbcmJRaLkX+BpWUvbmx6xf+ataXfHfy1WWvNxWuK5t++pEH2LN3N/cMeTBgGT9/P51njr+UZsWraPrs4/yStZdbTjyHNnmBa48AfV55mLcbd6d+yQaunPUZrZoewV3NO9J92yyabdrKC0dcQP7eeUw770a6bNhDsWRxy8wRjO49AoBnn3yQrVs3sfLwRrzf6AxoBDkfv8SHbQ+1/y7Lcu+EmtrhPLoUWuN/f658qLb80ZtT2bf/AFdc7z7warvU5oJPptJ57gLWHtc66N/0ZofO7BH3wc6zqrpfZuFrCN6kzX9yXY41FrmEDPa7jLRZY9864PXmZwfdv0pNmsxTRLhv07hhffmpa/g15FCNnDiU6W2O52yPzqKvmtdh+4j+VKpchUn5Z7FLqrF18igySorLEqinjX+th7p1WZ3ZnA/rb+DbaqeT8/M0xoeQzN9u3N0qI6MBy45qyaS6VpKbXusUqGWtU1Apj+cnT6D4KOs+G8+1+QeXFXzL08vm8NExVnLL33votNgE6AAD62q9ddVrl73u/erDNFy7medPuoi9Upk7fZxiF1TKo+CU4H8PwI6M6iGt52lGzY5cu906+yiWLJ8JX6WvDQdic82DJvM4GDP2Xp7oeA23z32XYf0Pdei8MfUpHndJ5L7GqgYy+sH7Kc7Joig7i1F3jPC5znvHdWJjRgM6V1/oNn9mjY7MPN09iZR2ro0muG/t4Wbzm7iPpR0z/j721KjCmNsOjYSY9OBQaH/oQpZAoxDW7XUfg/vo6t/43KVd17XZY3LeZUHjdB3V8G6j7tAo6CYh2ymRJXOAnjW9m5CUioYmc1vpSIbRvSO7eKJ0WJVr3fzcT1+i6sF91G9WH4D5hx1KfPdPGc3Uw917+pflHuoA6/vSQ+yqVIFWi1czeMihIWfDnxjOGce0p3O383my/VVl86d/8T4/nOneNhute58bS/Xtuxh6j/8hYusqWDEPnzScPVUq8nKHHgAc9syj9Ly1P4OeGcUL7d2vSJxRy//I1adPdF+3oMZRkYavVLmiydz2nN1jPxoYO7I/jRs144Z/3VW2fOyo/mBg0DDvjr9b3nocck8HQFzaw+ZWPB4qwmU7pgNQJJmMGXsvgwdN4Kfmh3mV49qu+2bp/SXqQsawuzCVcph08tVw7CUUbviSzriPJHAdQ9vl8zdZltWS12ULf4vVNut3VIYPrb/8jt1SFVqfB0DWqP4MHPoo+/Z5X4q8XWoy6uFB/Ked+9jlZ1q3oiewrGE9r20OivfwPH+2+rjIRCnlLaWS+cRxg9mZbQLWFKP14PC+TOrSk1N2FXAD8NKzj/Nwq6MoPPV6AAZ5rD922J1M63qT27xxw/rStMXh0NIab1uaRmdVzWdWx3x2Tx7FwjB68h/veiNXr/qs7PXmKj7uMORiUfaRAEzcuYZ91azRHC8cGfr+dnt06u2sazVqL121FJq4j5DYK5V5yiORA6zKagGEP6JCKRWZlLru/Mc2TXmq3RXWkLMIjXxsCPe84PtpJQDFmZkA/FrFSohr1qygMKO+13rT//cuQ58agfE4hAZ4vGtPHm7e2m2eq79y/T+FZfT4gTz7hPfd6V5zGd3wbbX2PPDIEK91Lv/I/b4Y31Vz7kLcseMG8Hina8LaZriPKwiVUrGRUjXz1ZWtU/aDPq6iC1Vpp1njsQPZnQVZe/dw//BDj4P6ue0xgPeVbp76ZVdjU5t/cFYD94cCvHC6Nc55Q8ahe0N81sj9arLNlf2X/WSHK0P4K+Dptt6df99W68C4kf2g8w0hlREqwfBSh3PD3u4/xzv/tBmllG8plcxLTc67DH/dlId/+TUVzAHmn3FmwDLmH96Qj+t2pqLZg+udojdVcK81i7jfOW3E48OYc1hLNtnt25+7XLoM3jevB9gp7vc3DnZVXzQedziRA3xyRD7bpZbj5SqlnJNSyXxHRuC2YrAS506BMaPvplGdhvS87R6f65XeEnWfy02GAJZluV8sIpLp9vqZcljbXJ/h4Hg+pVRMpFSbueu9Kzp98Z7X8odH3F02/USn6/i8ofUswrGj76HBzHmMemRw2fJZNUKrHa84ukWE0SqllLfG9j32nZZSNXNXy7JaMXbYnUzqehPNi1Zy8XczKMh3v0XMb9UO46KPX2TRKdZ47qfaHrpL8l6PGrkvDWbOA4+7zCmlVDRePC70W/GGI2WTOcC007oA1jC4SV1u9lq+JSOXLZV93CrOQyLunKaUKp+S9n7mibQ8K/SneyilVDKQGN1oK6WTuVJKKYtjyVxEMkVkrohMC762UkqVT7GplztbM+8LLAy6llJKlWNJncxFpAlwPvCcE+UppZQKj1M184nAAKDE3woi0ktECkSkYNMmZ57fp5RSyhJ1MheRC4BCY8ycQOsZY6YYY/KNMfm5ucGHCyqlVDpqHMW9pQJxombeCbhIRFYCbwDdROQVB8pVSqm0Uz0rM/hKEYg6mRtjBhpjmhhjWgBXAV8aY66NOjKllFIh03HmSimVBhy9nN8Y8xXwlZNlKqWUCk5r5koplQY0mSulVBrQZK6UUmlAk7lSSqUBTeZKKRUn1zasE7OyNZkrpVSc1MmJ3fOANJkrlaaOOrA40SGoONJkrlSaynC5713romUJjCS55O9Nz8dEajJXKk3llBSVTXdf8GsCI0k2sbqjeGJpMlcqTbnWzGP02EmVRDSZK6VUGtBkrlQUWhatSHQIKojeyz9ze20SFEesaTJXKgo1i3YmOoTQpEEGq1WyJaLtmmdXdDiS5KTJXCmVErTZPzBN5qrc2tA1L9EhxE1WUXHY2zQrXh2DSOJPkuBn4F9N6tKjYW3uaFYvZvtw4hmgFUXkZxH5VUT+EJGRTgSmlDqk3d7whxa2Xb6ybPrf194Z9vYN928Oe5tEqFeyMeBy49HGlIgWp1GHN+Gxo5pRLUaPjANnaub7gW7GmBOAPOAcETnZgXKVUrYrV2/noT+/Dnn9DV3zGN17eNnr3Pr1w95nsjWzi4+Irlrzhc/5QQqKSLKfyTnxDFBjjNllv8y2/yXb50CplFahYjbHtj0ppHVzzD6H9pr45glf+i78uGx64vX3hrTNsfsXxCqcpOFIm7mIZIrIPKAQ+MIY85MT5SoVK82LVjpSTti1Qj/+/dt7QdfJqVAppLKyORhtOEnKOtaVM92bKkJ5D6afc3XZ9GmF250NK0k4ksyNMcXGmDygCdBeRI71XEdEeolIgYgUbNq0yYndKhWx6sW7Ex2Cm+z9gRNw0cHwOzCj5S9FXlT4FXVKkr89vU6J7zwz8Ia+cY4kPhwdzWKM2Y71QOdzfCybYozJN8bk5+bmOrlblWaSeRRFuKfrZ237PkaRJE7bwj1c/Vvo7fdOC/XWBDlpe4bimxOjWXJFpKY9XQnoDiyKtlxVfmUa52uht899x/EygaC9Q4etWB+b/SZYnezKcd9nrFvw71s7izHb58d4L7HjRM28ITBTRH4DZmO1mU9zoFylonLc/j/KpptXrA1A4+K1QORt3fEaOdG0eE1kG4bgyjWfB13nqtX+18kAbr1jIJPW/cz5m79xMDJn9Vj+G02LV9O983khrV9dsqiRUyWsfWSbA5GEFhNOjGb5zRhzojHmeGPMscaYB5wITKlouSbeG27vz9Qdizh7Warcy9r9R6O4uMTtdTT3J3/8+gEBl9ctKWTiDYHXAbji2l48f3no49ePPPhnyOv60uCg1QZepWJoHcED/nU/s7tfRKMWh0e130B67Ps9ZmWHS68AVX6dsX1WokNw1DkXXxV1GcFq5rklhX6X3bLgQ695wS54SYRaxTsAMH4apyXC615y9/0d9jYbuubRd+E0Tt/xI8Or1eXfi6dx8/WRdWB22POL17wPMrZTt6yjNLXTYWpHr2KqyoHkOYWMhMTgaof8le53SRTgtnnv0nPxf+n7w2tcOv87t+WlIVQxO90u4imVE+Fpev1927zmXbPyf17zAv24xF2ETU4Dbx/CmxffymmdzmLkrUN8FBvaGz2140V8Xsu9U/Tk07sw6cABzt72Pddefh0lKXyJTOyeLqpS2u1z32Vv5RwI41YSNcx2/rb6wpNClo+O1CrbrFrn0YV/RVTm2FuHcuuKFfRYOpdlWa0AGN5vVNnykY96JxsIP495rn/9rX2Z/9vcstc37stha85iFuUcWTbv8F0lePrXr9+zqYJABFcvmjCjPn/zN3xct7Pf5VklgTu2K5vd7JHw2qwh9Dhr1a5Lrdp1veZ3O/d8upWVFV4yT6bLqrRmrnwa1n9U8JU8HLN7qd9lDUsiS57RqHFgj9e8wfc/xMCC97ir4yURl9usZUtOWRXagK1g48dD4euKzosuv56z5s31sba7O/uPZFTvEWHuMXBCy/STNnLtH0p/au7dG3D58m6dAocVgtIO7nhJpnq8JvNy4ugDi+jzw6uJDiMp9L33AVoebnWKXbZuOpWMd9KPVLX91nM3r1jzBYPnfkBFCXzye/LGxQDU8LHe0F/eB/zX/gbdN97tfiG5uQ0AOHqfMw9vdjJRNS9aFXB5v/kfRb2Pa1d9xiNRVpVLTHh/tSZzFXcVSw6Qsb8o+IouJMwHR5pkOucM0ZPX3sOxexeHv6F9bDzba/sPtBLspOvvpU//EcH3f809bOiaR597vc+Eala3mgRCTRj/vOom+s1+k3ty6oS4hW++mqdcSYbvtBFNH8V9fYZFvG3pe3BEcQ5dzrgg8iAAEySZ1zTefRXJQpN5OdH5p1+Sq4EvDmLRARpPWdnuw0aqVQ3ennzfgHF06X5hyPt4eOW3XqOWLl1oDyF0+POS4m8HAJVM4KaiRNJk7qL0AoBuf/+Q4EicVa9kIwNHTCx7Xc0EbtssFe6XL9WTZziyD1hnOVUPOHWHQrj9l3dc7lsen4N5bc8+VDro3q7f567AteTMzNDTRsfdc4DYnrU13m+N2KlZpXrI22SYYp/t622ahzcmPdCfVZetYZUVLU3mPhy7ML0e0uvUnf1i5agDgZs5Ll/7RYQlx+7vvvWca+mx6nMu3B5e01Ugw+4ezckLo7uwJhL+Em24o1k81+477z1OWhr+fXbCvYnXhKbH02/xJ1xx1c0hb7O8w9HMOu1Mr/l5x+dx2q7ZIZcT6BP2Y6cOIZfjBE3myq9kaAO/Y/abVN9lndo2Kl4X130Hui1tsyrPkIcAABseSURBVJYteezGAVxz4x2O7tN3ckiCNyICA/u5Xwxec2fwJor7F33C2D3hJfOjj2vLfbcOCmubilUqU6FCBZ/Lzjzge74vgZJ51TBvDRAtTeZhSMar9SIRbo0r5HIdLvb4/fMZMmBc2es6RYHvQ13F7HR0/yP7Ju7OFGIgMycnYft3SvO6jQDosP5P2mRXA6Ca8X8l6F23DeLiC6O/Ujcavf55PXOPTr07u2oydxGsOaLlvuA1w1N3/uxUOEnlrK3fec27cs3n9P3m/8peB7piVPyMErh2xadB9+3vsnJXNc02r5+oUNvwQyk/ct4X8gQWv1p4sOeKnrpgcXhXkPo43tdc+282dM3jiesOPRHIiWa/MWt/YNKW1L3DYSxoMnfxz9VfUb9kA5WDjA3257SdwR+wFOyihhpmO1ev+izkffq634Snixe4d+iG+mVyTXK1d3mPxc4uKmbg8MfKXufu8K4Z957zFrf+9q7ffdRZuzlo/vIXretdERd168pVC78MXFCUork3efgJzFAcw9spbOiax8fn3RBwncF3j+H3M87ymi9JkDZuvu42rrjs2kSHQf+TQrsjYzwk/l1JIrlrNvDrGedw1/1jI9q+UlFR0HpVsCaOLoWBr+rzHOdaY1/wC14e6DMy6DoQuE061CYUMe410aH3jGVE3/CvJvVZtkdC9Hzt694nTqp4sHw97MBTpnGqszf2Zx/RNokWl7h/juuXbPC5Xo0wRtDEmibzMPhrKihlCD5+osXeQw8rOHfLt2HVwgEuWuZey667c5efNcNXpSS8MbS+mieMRPCR8jhotYqsoZNNd24B/F+8VKE4cHINtTYc7H0tFY8xQTkHrL/psF3rk67NvCLWMMzMSN7jOPquHnx+5JHBVwzRndt+oqDzGT6XGWN4seFOvjgi9E7TWEnudyXJhFafCLzW0NrHOPYw4b5fv0idQu/OpA1d86hinEvy8dZ4xzb6zH6TPi3b+l3nwsKv6PSLM09cP+2XBW59HRXMPgdroeG5f+CD9P3pdfo1OCas7XLM/hhFdIhzP2ax/Vk87Jg8GjRp4Vh59154E9nZ2T6XlRjDuUedxnGNj3Zsf5GK+q6JItIUeAlogNXbM8UY83i05SaL2+a+w+QTL7NfhfchPPzgUpZkH+Y278ROHTnm7UmsqtsiwogOxTBwhHWYM0b15/FTr/e7Xvglh1aCY19JH79/g11GsTTadZCjDizm9Hm/81sH61nhz155FwATZ/p/2ESoJ/MDhk4AoIFd1iMbSzvW8kMswVkD738QgJUrQhtzfvOij6i+Yy90i8+45owsfzc0P/SJOGfLt0D4d2pMFpn2LQtyzD63RO55FmdMuB3cseNEzbwIuNsYczRwMtBbRNo4UG5SGN5/dNm06+iIYw4s9FrXag8/tNIF33/n+wHAUTQZ+tp04NBHQ1wzcRcP+RsxknfCSUG3vePukXx19pUMtpNcrF3W4xYu63GL1/xI3ra8tqdRzezgHysi7zwt7WepULGqz+VjbhvGffeN87nMWaEfgamX9YlhHLHn2WZeKtSRT+MbhzdW3glOPDZuvTHmF3t6J7AQaBxtuQkRZa7L8viVvm/kROrv9T+m1n+brv9Aqm11Yix1aH9ohstqPjtuHRjSd64DT/9xFU3H16m7ZgccbRTJx6PbWeexpFtnJtw8MOxts5P0cQMZSd5m7qRIP+E3HtG9bLp+TnzeR0ffFRFpAZwIeI3RE5FeIlIgIgWbNm3yXJzUOuz5he7bZoXUoVa5KMwRD2HecnPIfePDKz8MnrWOrIO+240TcW/yUN3w8wzalw7XDPPYvnPhv5jTPbq77qWm0FJWsJFYwcb1myB3Y0wFpc0spePvMzMCP0PvvbzD+Dzfuc7YQBxL5iJSFXgXuMsY7zs5GWOmGGPyjTH5ubnOXF1VpyQ+Pwofnn8Tr1x6u9tH2d/ntuo+q8e/cwhjzv2L/dCtNj6aiQAu2fAltwa4jB2g2T7fw7Ry10b2iDInr0i9e+BD1Nvt7JWg6aJeyUa6/P1j1OVkZkdX00zNmxNY2h2wRqO92Lge79TdT9Uqvpu+Sp1Sqyr1K/juPHWaI/V/EcnGSuSvGmMCZwIH/ePPH3j+qIvitTuOXLue744KsIJLTTDTbnOrvWsXBHiSWqw+2MHqpLf8tZ3lhW/wZHv3Zo7JPfqXTefsO3SWUcFnLd19L4OGPMKkAB2SoWq1xJl7sKRy0nCVjTPj238742xHygnGqihE1/nZZ8FHbK5WMepyImX8tJk/feYN3L9hDa2btopzRMFFXTMXaxDw88BCY4yvnri0Mea2IDfQ99GG3LNlp7je0+XCwq/8LnMdA351z9uQA4GH3w27+1Dn7+VHtiubDnVctqtA27g2X90+9x0GDpkQdvmuWqy0zhxar0ze5qBQXHz1jVyw6WtuLpiW6FAA6PHnDJoVr6Jdhy5u83st+NDtYdIZIVyPEczg3sN47PoBEUQZW9nZ2UmZyMGZmnkn4DrgdxEprZYNMsZ84kDZbqqZv9kpNQKuU8nsYa9UdnrXEWt76ikc/t9nKaxa3+fySBJjIKVD9gKJpMba/tSzwK51lz6Y11cbaWWziz0S+NQzECeOxpAB4xgC0DWPpx04U0ik567om+gQyoy9dSi+ro1+wL7ydvBk60pff31LIoHbl5OJv6cpJbOok7kx5jvidEbbrfAXPqzfNeA69YoLabh/Cz9WaRdwvZhwODEnq/ZzF1H5uIOcnFXLa9ktP31EYaO6vNHM+54eKjk5/an11weSDh2gySw5xz75EWot9ogNG/mxdQQ7iFEydvqWs7klhWzKqBdwnZ5zPuGp/CsCrlPaKlTR7GGfVOaUVaE9C3PAcPfWtN4Fb1k/513zGDTwIQDeCKNG7Hp80qWd2587f3ydgxWzmZx3WfCV01S6v8eJklLJPBRtNy+Py36a7N7KggTdOqPHt59wsEY16Oq/c2jovWN5yk9C9fzJqmp2sbLbKQHLC2TovZHdmMwJrYuWsSkzugcYx9OggdaFT5NTvPkn3fnrAE1maZfMb6h3Iu8vnxV8xSidta2Ew1a8w9Nt3WtYsb03tmWQy/M8o2FSsF3Q0/dn/jPg8hP2zefXisfGKZrUpDVl/5L9kYuuUv/b7OHkM7pyVFEFum+bFXS8tJcwEvG1N/dxG+3hVVSou0zAZ8XrIQ7xDyFiZ4Z5T/E792cyaPb7MYom/dSN4tqNVEp86SjtauYAPW+7h57AmAfvc7zsO757hQxTEnGTRDja75lLhz/+5ImTroz5vlLFy5f2Dmv98//RI0aRxF+DRo1gcWQXZgVSmoK7b5vFHQ1idxZTXFI+kn0DScwV7mmZzEvF4lmXQ4Y+7HiZ/giGwQPG8UTM2leDf7n6/PAaUnQQuuZxy4IP2V6zCql8NzzlX3ZJMSef0jnqcuLR1BhrVapUA7bRZv8SoGNI28zOb8B/18zjgqYnxDQ2f1Irmcf6Q5LAoYVHH1hExZIDtFiwkk2tGsV134FOjwcPeqhsOtZP8lHpLTMjdZJ87Vq1eaPyrxx/fOiPhWtarQG3tzknhlEFlnZt5onk76PafMsWt9e+zhiqFO/j03OvZ+AD/m8FPzjKtt8Mj4cL+3uCTyzUO7g1bvtyUqui+IyOUsmnS4eu1K5TN9FhhCy9k3moD66Msdw11uX8lYzH8zrDOBFoXrSSPgNCe5ZncPE/Azl77nxuXvQROWZfEkQTuvGZyXM1cSKFcqfM0s78DD8PbCgvbeaJklrNLEku2Ee19BmK4Rr9x6dUrxGDsdT2A4rrHdwSZMXo9Rtkja9++cto7iYZf527nVN2G4Py6staB6hcLfhVeP86/1qWzf2Mbut3B1xPR73ERlon84w9gT9UgVyy4Uvqbt0R+1ErIbTT33JH+A82CMWgYY+xe/ID1Ni4Dc65Oib7UCkgyAlsm7z2IRXTrHkrXmt+mwMBqUikdTIf9MAkCv77HLOqhv8sR9dbwYbK/aFx4Yl25M0tCz5kdYM6hDvSJOidIB2ntbJkpTXm1JbWyRwguyQxT1kP9WIgp1r1QxlpUpKkj/tKjp4NpVJbcn67/Yh5vSFJk53Tki15nvNTarWjp6tYXJfhVr52gMZU+checRLNaWr5+Zh7J4yBQx+hRdHK+IeiLAn88F234lNO3TU7cQGkEaceG/cCcAFQaIyJ212NJBl/6UMYu9162TqoA62W/8WO2tXjEJRSAcTpVE18XDQ04abYdO6XR07VzKcCibv0KSkd+qHxbD8fct94NnTNY8h94w/NTLa2jzhLwp/lckc7QFObIzVzY8w3ItLCibICSa185yNan8MQY/cFunHJxxRnZsTlpmCe7vjpNcjITMi+lSqP4jaaRUR6Ab0AmjVrFq/dxpchqR4dN77X4ITte8j9DwVfyUVq/VCrSGgHaGzFrQPUGDPFGJNvjMnPzc2N125Dlrfvd446woHm/jDvd9Jw4w7y9v3OaXPmR7/vFBDrERMqFWhSj4WUGmceyUcg1NTxv3Ovi6D0QEKLtt+QB+kHcK7Du08xqfD1Todbu6r0lVLJPNk1Xb+FjfVqJjqMlHDrb+9RlOWjPT8Vsnqa0bOl9OBIM4uIvA78ABwpImtF5GYnyk0lQ+d9wMg7nbqrYWwk4hF1/jQim9G9R3gvSOK8IknUHxITaf7npTtHkrkxpocxpqExJtsY08QY87wT5TohXp/P3v1GxHV/0UjGIWhJnMOVQ44/th0A3TbMTXAk6ancNLM0L1pJfuFS3m3UPe77NnaNrpKPe3mXN7XMNgqlPplZFdzmJ9/Pi3Ja+46nswGAPN4r57cVjoW0T+alNb78wmUJi2HQ8MdY9+ojNFm2Drqfn7A4ksFVP3zB1kZ1uLm3nyv/NKsrFZG0T+bJ4qlr7k50CElh0JAgD8RO4vaWZBjN0n3bLI5Yud7Ri7EabdoGdaDh5u2Olanir1wl8wyt9SWtxKfJ1PDKpbc7XuaY24bRYOID9Lkr3ve2V04qN3dNNALNl6+ny98/cvqOH2O+v2QaOVLKJGNQKiloIk996V8zd8lf9wx7BIAb330yQcEki+RL6skXkVKpJf1r5nr+nlo0qysVkbRK5r3mv+93mWsLQyzzheaiKJWTH99jDizkrG3fJzoMlUbSqpnlgT7JcwVmMl6Yk8xSIYc7eQXojLN7OFZWKqpTsoktGcl3w71UllY180BMKmSLGKtRxbpvzJnr5iQ4Em+p8NOXDEMT08UzGcU8uUkvHHJSitfMDa2KlrM8q1WgVUKZVS70vns4vQFI4gdGlNc3p5w5ras+mMxpKV8zP//rLxMdgpuau/YCUGvXngRHopQqT1I8mQuDRyfXMMMrjzydXvPf54aWpyQ6lNSkLRlKRSTFm1mST7vTOtHutE6JDkMpVc6keM08QtqRlby0zVypiJTPZK6UUmnGqScNnSMii0VkqYjc70SZvnhX2rQal3b0pEmpiESdzEUkE3gK65HEbYAeItIm2nKVUkqFzomaeXtgqTFmuTHmAPAGcLED5arySE+2lIqIE6NZGgNrXF6vBTp4riQivYBeAM2aNXNgt6EJdll9z8X/ZXuNKiT1hTRKKRWEE8ncVyunVwY1xkwBpgDk5+fHrf5lgjTCjrt1aJwiUUqp2HGimWUt0NTldRPgLwfKVeWRdoAqFREnkvls4HARaSkiOcBVwEcOlFtu5S5fT+ui5bT/dVGiQ1FKpYiom1mMMUUicgfwGZAJvGCM+SPqyHzxuNgnlEpcKt6KduADExkIcOaliQ4l/lLv7VIqKThyOb8x5hPgEyfKUkopFT69AlQppdKAJnOllEoDqZXMHXxsl1LhqnDwYKJDUMqv1ErmHrL3FUW0nf4kqHB02jUbgKzikgRHopR/KXk/8y5//0jr9YUMu3t0okNR5UD1/fugKlTcrzVzlbxSK5nbQxOzS4oZc9uwBAejnFS7aAersiDjYHGiQ/FydeUm1FnxKbef0SPRoSjlV2olc1skFwm6Plk92CX+Kv7OnPs7x9Vby5D7xyU6FC9nnncJZyY6CKWCSMlkHg5N26mh/6AJiQ5BqZSW0h2godDOTqVUeZD2ybyU6LBGpVQaKzfJ3JXRthelVJopN8nciGZwpVT6SqlkHklDiaZwpVR5kFLJXCmllG+azJVSKg1ElcxF5HIR+UNESkQk36mgnFSxyLoEO+dgZPdxUUqpVBDtRUPzgUuB/zgQS8TO3/w166rWBvK8lh3353oqtPySE/dmxz8wpZSKk6iSuTFmIYAkeKTI85f39bus36DkuzxcKaWcFrc2cxHpJSIFIlKwadOmeO3Wpxp79wLQvGhVQuNQSimnBE3mIjJdROb7+HdxODsyxkwxxuQbY/Jzc3MjCrb5n2vJMfs58s/VEW1fqvRq0BZ71kdVjlJKJYugzSzGmO7xCCQUQ4Y+zBCAbh2iKqde4XZoBC03bnEkLqWUSrS0v2uiL8P7jSZn6J0MHDUp0aEopZQjoh2aeImIrAU6Ah+LyGfOhBV7msiVUukk2tEs7wPvOxSLUkqpCOkVoEoplQY0mSulVBrQZK6UUmlAk7lSSqUBTeZKKZUGNJkrpVQa0GSulFJpQJO5UkqlAU3mSimVBjSZK6VUGtBkrpRSaUCTuVJKpQFN5koplQY0mSulVBrQZK6UUmkg2odTTBCRRSLym4i8LyI1nQpMKaVU6KKtmX8BHGuMOR74ExgYfUhKKaXCFVUyN8Z8bowpsl/+CDSJPiSllFLhcrLN/CbgUwfLU0opFaKgzwAVkelAAx+LBhtjPrTXGQwUAa8GKKcX0AugWbNmEQWrlFLKt6DJ3BjTPdByEbkBuAA4wxhjApQzBZgCkJ+f73c9pZRS4QuazAMRkXOA+4DTjTF7nAlJKaVUuKJtM38SqAZ8ISLzROQZB2JSSikVpqhq5saYw5wKRCmlVOT0ClCllEoDmsyVUioNaDJXSqk0oMlcKaXSgCZzpZRKA5rMlVIqDWgyV0qpNKDJXCml0oAmc6WUSgOazJVSKg1oMldKqTQQ1b1ZFNz05zSMMdA1L9GhKKXKMU3mURr77yGJDkEppbSZRSml0oEmc6WUSgOazJVSKg1ElcxFZJSI/GY/ZehzEWnkVGBKKaVCF23NfIIx5nhjTB4wDRjmQExKKaXCFFUyN8bscHlZBTDRhaOUUioSUQ9NFJExwPXA30DXAOv1AnoBNGvWLNrdKqWUciHGBK5Mi8h0oIGPRYONMR+6rDcQqGiMGR5sp/n5+aagoCDcWJVSqlwTkTnGmHyfy4Il8zB20hz42BhzbAjrbgJWRbirusDmCLdNhFSKN5VihdSKN5VihdSKN5VihejibW6MyfW1IKpmFhE53BizxH55EbAolO38BRPiPgv8/TIlo1SKN5VihdSKN5VihdSKN5VihdjFG22b+XgRORIowapp3xp9SEoppcIVVTI3xvzTqUCUUkpFLhWvAJ2S6ADClErxplKskFrxplKskFrxplKsEKN4HesAVUoplTipWDNXSinlQZO5UkqlgZRK5iJyjogsFpGlInJ/gmJoKiIzRWShiPwhIn3t+bVF5AsRWWL/X8ueLyIyyY75NxFp61LWDfb6S0TkhhjGnCkic0Vkmv26pYj8ZO/3TRHJsedXsF8vtZe3cCljoD1/sYicHcNYa4rIOyKyyD7GHZP12IpIP/szMF9EXheRisl0bEXkBREpFJH5LvMcO5Yi0k5Efre3mSQi4nCsE+zPwW8i8r6I1HRZ5vOY+csR/t4XJ+N1WXaPiBgRqWu/js+xNcakxD8gE1gGtAJygF+BNgmIoyHQ1p6uBvwJtAEeAu63598PPGhPnwd8CghwMvCTPb82sNz+v5Y9XStGMfcHXgOm2a/fAq6yp58BbrOnbweesaevAt60p9vYx7sC0NJ+HzJjFOv/AbfY0zlAzWQ8tkBjYAVQyeWY3phMxxboDLQF5rvMc+xYAj8DHe1tPgXOdTjWs4Ase/pBl1h9HjMC5Ah/74uT8drzmwKfYQ3VrhvPY+v4lzFW/+w/7DOX1wOBgUkQ14fAmcBioKE9ryGw2J7+D9DDZf3F9vIewH9c5rut52B8TYAZQDesO1sK1tVnpV+SsuNqfwg72tNZ9nrieaxd13M41upYCVI85ifdscVK5mvsL2KWfWzPTrZjC7TAPUE6ciztZYtc5rut50SsHssuAV61p30eM/zkiECfeafjBd4BTgBWciiZx+XYplIzS+mXp9Rae17C2KfKJwI/AfWNMesB7P/r2av5iztef89EYADWhV0AdYDtxpgiH/sti8le/re9frxibQVsAl4Uq1noORGpQhIeW2PMOuBhYDWwHutYzSF5j20pp45lY3vac36s3IRVQyVITL7mB/rMO0ZELgLWGWN+9VgUl2ObSsncV5tRwsZVikhV4F3gLuN+K2CvVX3MMwHmO0ZELgAKjTFzQogn0LJ4HfssrFPXycaYE4HdWE0B/iTy2NYCLsY6zW+EdQvocwPsN9HHNphw44tb3CIyGCgCXi2dFWZM8fg8VAYG4/uZDnGJN5WS+Vqs9qhSTYC/EhGIiGRjJfJXjTHv2bM3ikhDe3lDoNCe7y/uePw9nYCLRGQl8AZWU8tEoKaIlF7967rfspjs5TWArXGKtXT/a40xP9mv38FK7sl4bLsDK4wxm4wxB4H3gFNI3mNbyqljudae9pzvKLtT8ALgGmO3OUQQ62b8vy9OaY31w/6r/X1rAvwiIg0iiDeyY+tU21ys/2HV2pbbB6y0c+OYBMQhwEvARI/5E3DvWHrInj4f986Pn+35tbHah2vZ/1YAtWMYdxcOdYC+jXtn0O32dG/cO+nesqePwb3DaTmx6wD9FjjSnh5hH9ekO7ZAB+APoLK9//8D+iTbscW7zdyxYwnMttct7aQ7z+FYzwEWALke6/k8ZgTIEf7eFyfj9Vi2kkNt5nE5tjFJHLH6h9Ur/CdWj/XgBMVwKtYpz2/APPvfeVjtcjOAJfb/pW+KAE/ZMf8O5LuUdROw1P7XM8Zxd+FQMm+F1Vu+1P6QV7DnV7RfL7WXt3LZfrD9NywmilELIcSZBxTYx/cD+0OelMcWGIl1p9D5wMt2ckmaYwu8jtWefxCrtnezk8cSyLf/9mXAk3h0XDsQ61KsNuXS79kzwY4ZfnKEv/fFyXg9lq/kUDKPy7HVy/mVUioNpFKbuVJKKT80mSulVBrQZK6UUmlAk7lSSqUBTeZKKZUGNJkrpVQa0GSulFJp4P8BE/Yxw/J5s8IAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "continue_train = True  # set to True to continue to train\n",
    "num_big_steps = 70     # number of small steps\n",
    "num_small_steps = 200  # number of big steps\n",
    "if continue_train:\n",
    "    for k in range(num_big_steps):\n",
    "        for j in range(num_small_steps):\n",
    "            mine_dXY_list = np.append(mine_dXY_list, np.zeros((rep, 1)), axis=1)\n",
    "            mi_list = np.append(mi_list, np.zeros((rep, 1)), axis=1)\n",
    "            for i in range(rep):\n",
    "                minee_mine_list[i].step_mine()\n",
    "                mine_dXY_list[i,-1] = minee_mine_list[i].forward_mine()\n",
    "                mi_list[i,-1] = mine_dXY_list[i, -1].copy()\n",
    "        # To show intermediate works\n",
    "        for i in range(rep):\n",
    "            plt.plot(mine_dXY_list[i, :],label='dXY')\n",
    "            plt.title('Plots of divergence estimates')\n",
    "        display.clear_output(wait=True)\n",
    "        display.display(plt.gcf())\n",
    "    display.clear_output()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save current results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Current results saved.\n"
     ]
    }
   ],
   "source": [
    "overwrite = False  # set to True to overwrite previously stored results\n",
    "if overwrite or not os.path.exists(chkpt_name):\n",
    "    minee_mine_state_list = [minee_mine_list[i].state_dict() for i in range(rep)]\n",
    "    torch.save({\n",
    "        'dXY_list': dXY_list,\n",
    "        'dX_list': dX_list,\n",
    "        'dY_list': dY_list,\n",
    "        'mine_dXY_list': mine_dXY_list,\n",
    "        'mi_list': mi_list,\n",
    "        'minee_mine_state_list': minee_mine_state_list\n",
    "    }, chkpt_name)\n",
    "    print('Current results saved.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate the ground truth mutual information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ground truth is 4.982193620464953 nats.\n"
     ]
    }
   ],
   "source": [
    "mi = g.ground_truth * d\n",
    "print('Ground truth is {} nats.'.format(mi))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Apply moving average to smooth out the mutual information estimate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "mi_ma_rate = 0.01            # rate of moving average\n",
    "smooth_mi_list = mi_list\n",
    "for i in range(1,smooth_mi_list.shape[1]):\n",
    "    smooth_mi_list[:,i] = (1-mi_ma_rate) * smooth_mi_list[:,i-1] + mi_ma_rate * smooth_mi_list[:,i]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the mutual information estimate after different number of iterations. The red dashed line shows the ground truth, and the green dotted line is the number of iterations where 90% of the ground truth is reached."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3gUVffA8e9NCCkECBB6iEGUIi1CKKEZpEpTEQFFBAv4E7sGK9i7+ArYXnktICKggCLShEhQIAgBQydICRJ6Qg1JIOX+/tjNsptsGuzubLLn8zz7MDszO3N2suzZuXfmXKW1RgghhMjPy+gAhBBCuCdJEEIIIeySBCGEEMIuSRBCCCHskgQhhBDCrgpGB2AtODhYh4WFGR2GEEKUGZs2bUrRWtd0xrbdKkGEhYURHx9vdBgCOHToEAANGjQwOBIBcOis+e9RVf4ewpZS6qCzti1NTMKumJgYYmJijA5DmI38aSQjfxppdBjCw7jVGYQQwr4J3SYYHYLwQJIghCgDel7b0+gQhAeSJiYhyoD9p/ez//R+o8MQHkbOIIQoA+5feD8AsaNjjQ1EeBRJEEKUAa9FvWZ0CMIDSYIQogy4Kewmo0MQHkj6IIQoAxJTEklMSTQ6DOFh5AxCiDLgoV8fAqQPQriWJAhhV48ePYwOQVh5u8fbRocgPJAkCGGXlNhwL50adDI6BOGBpA9C2HXo0CFLPSZhvO0ntrP9xHajwxAeRhKEsEtqMbmXR5c8yqNLHjU6DOFhpIlJ2DVgwACjQxBWPuj1gdEhCA/k1AShlEoCzgM5QLbWOsKZ+xOOExwcbHQIwkq7+u2MDkF4IFecQXTXWqe4YD/CgRITTdfcN2nSxOBIBEDCsQQAwuuEGxyJ8CTSxCTsiouLAyRBuIsnlz0JyH0QwrWcnSA08JtSSgNfaK2n5V9BKTUWGAvQytcXoqJsVxg6FMaNg/R06Nev4B5GjzY9UlJgyJCCyx9+GIYNg0OHYKSdAVeeeQYGDoTERHjooYLLJ0yAnj0hIQGefLLg8rffhk6dYN06ePHFgssnT4bwcFi5Et58s+DyL76AJk1g0SL48MOCy2fOhAYNYO5c+PzzgsvnzYPgYJg+3fTIb8kSCAiAzz6DH34ouDw21vTvpEnw66+X54eHg7f35edvvAH5O61r1ID5803TL7wA5qRiERIC331nmn7ySdMxtNa4MUwzfyTGjoU9e2yXh4ebjh/APfdAcrLt8shIeOcd0/Qdd0Bqqu3yHj1g4kTT9C23QEaG7fIBAyA62jSd/3MHbvXZmzwvzTRvulWc5fWzB+DvD0uXmqbls1dwufVnz4mcnSA6a62PKKVqASuUUru11n9Yr2BOGtMAIipX1k6OR4gyKTwt0OgQhAdSWrvmO1kp9SqQprWeVNg6ERERWsakdg/Tzb8IR48ebWgcwmTj4Y2AdFaLgpRSm5x1AZDT7oNQSlVSSlXOmwZ6A3KnjxBXYPyK8YxfMd7oMISHcWYTU23gJ6VU3n6+11ovc+L+hCi3Pun3idEhCA/ktAShtd4PtHbW9oXwJC1qtTA6BOGBpNSGEGXAukPrWHdondFhCA8j90EIUQa8GGO6jFXugxCuJAlC2CW1mNzLFwO+MDoE4YEkQQi7pBaTe2kSLHe0C9eTPghhV2JioqUekzDe6qTVrE5abXQYwsPIGYSwS2oxuZdXYl8BpA9CuJYkCGHX0KFDjQ5BWPn61q+NDkF4IEkQwq6AgACjQxBWrq12rdEhCA8kfRDCroSEBBLyV8AUhlm5fyUr9680OgzhYeQMQtiVlxzCw2WAGnfw5h+mct09r+1pcCTCk0iCEKIMmHn7TKNDEB5IEoQQZUCDqg2MDkF4IOmDEKIMWLZ3Gcv2SjFk4VpyBiFEGfDumncB6HtdX4MjEZ5EEoQQZcCcIXOMDkF4IEkQQpQBdQLrGB2C8EDSByFEGbAocRGLEhcZHYbwMHIGIUQZ8GHchwAMbDLQ4EiEJ5EEIeySWkzuZd7QeUaHIDyQJAhhl9Rici/BATI+h3A96YMQdkktJveyYNcCFuxaYHQYws1orZ26fTmDEHZJLSb3MvWvqQAMbjbY4EiEO1m1/YhTty8JQtg1evRoo0MQVhYOX2h0CMKFlv39L03rVyOsVuUi1/vrnxNOjUMShBBlQFW/qkaHIOzIvJTN5MXbeGpAK3x9vC3zT5zNoGIFL+L3neSDhVsY1rkRo6Ka4O2lCt3WOwv+JnbH5TMCH28vfn3xlkLXz8rJtVnfGSRBCLvWrVsHQKdOnQyORADM3T4XgGEthhkcicgzY1Ui36/ZC5iaepZP7A/A58t38POGJJt1567dx/mMLJ7o3xKtNbla89a8zXRuWoesnFxublm/wJd9Vk6u3f2ez8ji923JfLZ8p+PfVD6SIIRde/bsASRBuIvP4z8HPCdBzIvbT8fGtQipEWh0KBbHzqSzeNO/RDapzem0i5bkkGf9nuOkX8wukBzyLNn8L0s2/2szb23icQA++nWb/X2eTudsxiWur1sVL2U6+xgy6berfCclJwlCiDJgyYglRofgdLlac8ubl9/n/1busvwqL6msnFwGvL0UgAd6NGVop0YF1pm/fj/bDp5i84EUZj3Rg8r+PkVuMzsnl+UJh5i6ZDsAP6zbZ3e9V+bG253/9MBW/GfR1hLFX8m3Ague7cN9n67iyKl0Rn2yymb5R/fZ/mD78uGbCH25RJu+IpIghCgDAnzK/30pw/9TcEjVF7/fwPO3h5NyLpM6QQEE+Bb+lXUpO4eB71wuif5VzG6+itkNwKioxtzd9XoWbzrItBW7LOsMmfRbgSTU543FACyf2J+Ne08wYfbGUr+XMT2bMSTy8jjiSSfOs+CvA8W+bsZjNwPQNzyUr3/fXWD5U9+ss3neINi5Z1iSIIQoA77b+h0A97S6x+X73pl8mo3/nODeqMYoVXgna57U85lU9vehYgVvm/nbDqYS/e16AipWIP1SNssn9udiVg67D5+hWUgQZ9MvFdjWpn0nuXPSCsvz/4yOpHmD6nb3++O6/YXGNCN2D9UCfS1nAYU5fibdMp2dk1tsclg+sb8loeSZ/VQPqgf62cx7qPcN9G4dQkhwID7eXvy+7TAhNSoRXMWPuz6KAWDZhH6W43tLmwZ2E0Ser8dFUb9GpSJjcwRJEEKUAV9u/hJwTIK4mJXDoHdNv7R/eKYXVQMqFrru1zG7mWtuUsnRmlFRjTl8Kp3QQn65pmVmcffkGMvzz8d25draVej/1hKyc003daVfygZMX8b3fryqwDZ+jO7Fyi3JfGH1Sz/P09PjLL/4F8UfZPGmg/z3oW4cPZ3Ot6v3FPm+J+dr5/ev6I2vjzdaa8sX8/iZ6y3LP15qm0zmPt2TSn4+pJ7PZNTHq7itfRgAH9zbkfHfrqd2VX8+uq9TgeSQp2HtKpbpm1vWt0y/P7Ijl7JzbJJvFf/Lf5Pqgb4M7dSI//52uVPaFckBQDn7TrzSiIiI0PHx9tvxhGtNnz4dkPsh3EVWThYAPt5Ft5fndyotk7s+iqF+9Up8OCqS7YdO8ea8zTbr5G9i0VoTs+0wHyzcUuh2JwxpQ9dmdW3m5Z0h5Ld0Qj+bvoWifDamK43qmL5ID6de4P7PYguskxdv3i/314dHsDwhmbW7jwHw/ZM9mLt2Hws3JgFwZ+S1/Bh3+eyiQY1KfPJgF77/c68l+T3RvyU9W9W3aaLKE1azMh8/2LnAGZGrZV7K5tb3lvPFQ91s7o9QSm3SWkc4Y59OP4NQSnkD8cBhrfUAZ+9PiPKoJIkh41I2R06lUyfIn8Ef2F7pcvjUBYZ/VLCNHyAnV7Pv2Fke+2ott3doyE8laCvPSzI9WtanSb2qRV5yWdLk8NSAlpbkAKZfydbJ69mZ69nx7ykAfrTqKH55zuUflT8924cA3wrcd3MTalbxo8+NDajiX9EmQUx7+Ca8lGLltmTLvCmLt3EoJc1uXFMeMD45APhVrFDqTvur5YompieAXUCV4lYUQtj3341fsTP5NB8OeprU85nUrOLHJ0t3MLp7E0sT0W3vLb+ibb81f7Pl17e95PDasAhah9Xgq5jdLIo/aLMsZtthYrYdtpl3c4t6bNh7grTM7ALbemdEB8JqBVra3fN35hZlS1IqQIE2f2t5ndj+FStwp9UVTD8924e//jlOVPN6lqacHi1DbK5IyutE/s/oSJ6eHmeZ7+djfHIwilMThFIqBOgPvAU87cx9CVGevbhkMgBZJzuw//g5y/wlm//l1xdvsVzaaW1A21B+3fQvQyKvZZ7VL+inB7aiZ6v6fPN7Ij/G7bckh/zCalbmi//rZnn+6C0tePSWFvy84QCfF3LG8MrQtnRqYhr9Lv8X+acPduG6uqY7wmc/1YPk1Au0uqZGSd5+icwf37vQZQG+Fejeor7NvAd6NOWOjg0Zlu/qqeYNqhM9qDWTftnCgmcL36YncPYZxGTgWaDQgiJKqbHAWIDQ0FAnhyNKSvoeXOeR//3J3mPn+Pm5PvhXvPxf0vpyy3Ze7wLYJIc8q+2UW2hYqzKP9WvJY/1aknEp25IgrJsorqtTsHxHrar+nDibYWmqsee29g0Z1C6MlVuT+fAX2+v7b2x4uSz5W3e356XvNxDVvB4vDL7RZr3qgX6FduYWJv8VQ4te6MsfO4+yce/JAtsvqaBKvnbn92odQq/WIVe0zfLEaZ3USqkBQD+t9TilVBQQXVwfhHRSi7Jg076T1Krqf1XXoJ/PyOJ/K3eyPCHZ7vJHb2nBJ+araPJ+zZZEoF8F5o/vU2D+F7/tpHVYDTo2rm2Zt+/YOcb970/L8/dHdqR1mON+0YOpX8Q66TnCoZQ0agf5O6xfwPo4uLqN3xHKaid1Z2CQUqof4AdUUUp9p7V2/YXcotSkFpN9WTm5vPj9BuDKv0w270/hhVl/FbnOt7GJlulJv2whWZuurglRfQG4pmYgI7s15s35l69I+nxs10Krfz7U+4YC8xoE214qeU1Nx9905ejkAI6/OaxRnSqMv7U1kVbJU5g4LUForV8AXgCwOoOQ5FBGJCfb/2XrSfYfP8d3f/zD7e3DmLx4G68Ni2DPkbOFrp+Tq/FSoJRCa83Bk2nUrOpHJd/LVyAdSkkrNjkAnMvIsnl+TP+Jn483ITl9LU1R2VbF3KY+0Jlra5fuOpD8v8CLuh+ivOvZSpqT7HHJfRDSxCTKoqKulsmTdyNYdk4u/c0dxc/dFs7kxdu4mJUD2N5ZO/TDFTZ3DAf4VuCnZ/vwb0oaYz5fXeh+WoRW58NRkVfzduyyfo9lsXlFlN0mJgutdSwQ64p9CZHfufRL+FX0Zsribazcarokc+oDnWlSL8ju+tYF34rz8LQ/WT6xP1sOplrmvfez7VCtd30Uw/KJ/cnV2iY5TL6vE81CqgEQGhxo+YLOu8P4rbvbs37Pcc6lXyL61tYlf8NXYFC7a5y6fVE2SakNYdfKlaZL/3r27GlwJFfvzg9XFJi3ZNO/NKkXxOm0i1QJ8MHbyzQ8+/mMrFKXU06/mM2LszYUuU7fNxZjfa5e1K/12kEBluURjWoC8NnGzwAY125cqWIrzoM9m/Llyt2M69PcodsV5YMkCGFXWeiDOHYmnVEfryKkRiW+GhdVYHn+8tHWliUcYmDENTzy5RoAXh0awas/FGzeXDqhH+t2H+PGa4NJy8hCKUWtqv6WEhYAt79v/wa1CXe0sXQiWyeHqOb1SvEuTRbtWQQ4PkHcGdmIOyMLlsQWAiRBiDIqOyeXUeZCb8mpF0g9n0n0t3EcOZXO8M6NuO/mpkTPiCvwurlP97TcGJWXHIACyeHWdmGM62v6Vd3FXHPIurO5eqAffcMbsCzhkM3rlk7ox4mzGdSu6o9SCu+fFDm5l9PDf0ZHWpqVSmPpiJI1eQnhSF5GByBEUX7fdpg+byymzxuLOXE2wzI/f0mIuyfHcOSUqVTznLX70Fqz49BpAGpU9uXH6F78/Fwfgir5ElK9+EqYecmhKE8NbGWZbhYSxIzHuuOlFHWCAizlHAa0vdy2//nYrjRvUN0yMpgQ7k6quQq73KGa6ytzNrL+nxOW50GVKjL36V7A5atvKngpSxnpwuRv7z948jxj//uH3XW7NK3D4/1blviSz1XbD1Otki/hVncQW8vOySX1fCa1g65uwJ8p66cA8ETHJ65qO6L8ceZVTHIGIQyXdOK8zYDtWmvW7j5mkxwAzlwwXQGUkJRimffh6NLfyHdNzcs3ky2b0M8y/dxt4Uy8s22p7gfo3qJ+ockBoIK311UnB4CYAzHEHIgpfkUhHEj6IIRL7Dt2lsQjZ+nXpmC9rYe+MP2ab9GgOsFV/Bg59XdOnsu0u50vVuxkwXpT89I7IzrQtH4QXZvVZUjktQQFVKRaoK9lMBwwNevYs3RCP3JyTQPFtLqmOlsPnqLD9bWu9m06zS93/WJ0CMIDSYIQBVzMykFrcFRT+brdx3jtx02W5+2vq8X/Vu4ibs9xRna73jJ/xJQYFozvbZMcXh0aQWST2jw9fR07Dp22JAeAFqGmzt4JQ9rY7O+X5/sy6N1l3N31ukLvLvZSCi9v0xt8bVg70jKzqORXusF4hCjvJEGIAga9u4zmualcX9f+jWSlZZ0cpiy2HfbxyxjbcXcPpdoO2hLZxFQfJ6/D2Vphxdp8fbxLdVdwgG+FQiuXuotJ6yYBEN0p2uBIhCeRPghhVzYV2HXU/ghbpXHmwsVSrf/E1+vszn9/ZEfL9H3dmzDnqbJ/A19pxCXHEZdc8LJdIZzJvX82CZdJOnGeh774gxmPdgcg0es6AI6cukC9ElwWWpi8kblCgwP5t5AhHb98+CZWbEm2jA8MpjGCrZuHWofV8OhaQfOHzjc6BOGB5AxCcCk7x9JRPOqTVTbL7vs0tlTbyrtnIc/hUxcA+ODejvh4X/64je3VjOmPdue7J26mQXAgo29uYrOdfm1CaVrfMU1cQogrIwlC8PC0PwvMC81NJjTXVG7D3r0y5zOyuO29ZZyzKj5nvV6fNxZzOu1y81JQJV8WvWAay+ChXs24o+O11K0WQM0q/oCp0/iVO9sCuH1/gBHeXfMu76551+gwhIeR/4kebv2e4ySnXigwv7q/F+cyTF/+K7Ymk3Iuk2Yh1bixYTDZObmWgnZ3frjCMhLZ7DV7bbbx8pyNNs+VUkU2E3VqWodvHom6qiat8irhWELxKwnhYHIG4eFemWv/zvUnHryb/V5hAGw7eIoZsXt4/jvTQDf985XCfnbmeg6nXmBG7B6b+XuOmgbX6ViKkbokOdg3Z8gc5gyZY3QYwsNIgvBg1k1C86J709s8SPvQTo2oV70SPz9nGtv4ty2XK7vmFFLW4v7PYgvdz3gnj2UghHAOaWLyYH/uOmaZruzvwzODWvPMINOX+aJFi8xLbH9DPGI1yH2dIH+OncmgOIFyA9pVe2P1GwBMvGmiwZEIT1LsGYRSqrFSKkYptd38vJVSaoLzQxPOdCotk7fMYxXYG8oyNTWV1NRU7oy81mb+gRPnAejcpDYzHruZ+eN7F3jtRKs7m+0tF6WXmJpIYmqi0WEID1OSM4j/AeOBLwC01luVUt8DbzozMOFceYPdAITWDCx0vf3Hz9mdP/KmxkDBs4MlL92Ct5cXnzzYhY17T8jZg4N8N/g7o0MQHqgkCSJAa71B2RbmyXZSPMIFUvIVwqtcxJf4s7eF88Bnq/nm0SjunHR56M6GVjexLXnpFv45epYalf0sQ3deX7cq19et6uDIhRCuVJJO6hSlVCPMoyYqpYYAR50alXCa42fSGTHl8tnD8on9UUVU5Quq5Mv88b2p4l+Rl+4wNR2Nimpss463lxdN61ez3NMgHO/lVS/z8qqXjQ5DeJiSnEE8AkwDmiqlDgMHgBFOjUo4zb0fX75Teu7Tpatn1O2GunS7wXPLXRjp0LlDxa8khIOVJEForXVPpVQlwEtrfV4p1dDZgQnH++9vO22eB1XyNSgSUVrf3PqN0SEID1SSJqb5AFrrC1rr8+Z585wXknAW63GcJ93bsYg1hRCiiDMIpVRToDlQVSk12GpRFcDP2YEJ5wmoWIFmIdWMDkOUwgsrXwDgnZ7vGByJ8CRFNTE1AQYAQcBAq/nngTHODEo41sqtyXywcIvl+U/mO6SLUqNGDWeGJEopNSPV6BCEByo0QWitFwILlVKRWmsZqaQMs04ODWqUrNbRwIEDi19JuMy0gdOMDkF4oJJ0Uv+tlHoEU3OTpWlJa32/06ISDvPN77ZDen42tqtBkQghypqSdFLPBOoAfYDVQAimZibh5g6nXmDO2n028wobxzm/RYsWWdVjEkaL/i2a6N9kPGrhWiU5g7hOa32nUupWrfUMc5mN5c4OTFw96wqrrw6N4Lq6VQpfOR9/f7npzZ1kZBVfFFEIRytJgsgy/3tGKdUCOAaEOS0icdUyL2XbjM1wd5friGxS8jEZAHr2LN1NdMK5Pu3/qdEhCA9UkgQxTSlVDZgI/AIEAsXe86+U8gP+AHzN+5mntX7lKmIVJZCTq5keu8fmnod7bmpcxCuEEMK+YhOE1vpL8+Rq4Nqi1s3nInCz1jpNKeUDrFFKLdVar7+COEUJvP7jJtbuPmYzb+7TPfH2KrzWUmF++OEHAIYOHeqQ2MTVeXLZkwBM7jvZ4EiEJyk2QSilgoB7MTUrWdbXWj9e1Ou0abiyNPNTH/PD/nBkwiHyJ4elE/rhVUQhvqKkp6c7IiQhRBlWkiamJcB6YBuQW5qNK6W8gU3AdcCnWuu/7KwzFhgLEBoaWprNCyu5VsOH+vl4M2lU5BUnB+F+5MxBGKEkCcJPa/30lWxca50DhJvPQn5SSrXQWm/Pt840TNViiYiIkDOMK3TLm0ss0wuf72tgJEKI8qJE90EopcYopeoqparnPUqzE631GSAWkG8uJ8jJvZxXo5rXMzAS4SyPLH6ERxY/YnQYwsOUJEFcAj4A4jA1F20C4ot7kVKqpvnMAaWUP9AT2F30q8SVsL5i6dnbWhsYiXAWfx9//H3k3hThWiVpYnoa081yKaXcdl1ghrkfwgv4QWv9a2kDFMX7/s9/AHh/ZEfLkJ+ifJnUe5LRIQgPVJIEsQMo9SUtWuutwI2ljkiU2oWLpiHCW4dJBVYhhOOUJEHkAAlKqVWY7m0Air/MVTjf8oRD/GfRVqPDEC4wdtFYQKq6CtcqSYL42fwQbmT++v1MW7HL8rzD9bUcuv2QkBCHbk9cnRr+cnYoXE9p7T5XlkZEROj4+GL7vz1ertY2l7XC1d0UJ4Qou5RSm7TWEc7YdlFDjv6gtR6qlNqGnTugtdatnBGQKFr6xWxuf/9yMd3b2ofxcJ/mBkYkhCivimpiesL87wBXBCJKxjo5AE5LDlKLyb3ct/A+AL659RuDIxGepKghR4+aJ8dprZ+zXqaUeg94ruCrhDP99c9xm+fLJ/Z32r6kD8K9NKjSwOgQhAcqtg9CKbVZa90m37ytzmhikj6IovV5YzEAI7tdLyW8hRCAcX0QDwPjgEZKKetrKSsDa50RjChcyrlMy7QkByGEKxTVB/E9sBR4B3jeav55rfUpp0YlChgxJcal+5s+fToAo0ePdul+hX33LLgHgO8Gf2dwJMKTFNUHcRY4q5SaABzTWl9USkUBrZRS35oL8AkXWzahn9EhCAM0qdHE6BCEByrJjXLzgQil1HXAV5iGHf0ekG8qF0g9n8ndky+fPSi518EjTbxpotEhCA9UkspuuVrrbGAwMFlr/RSmQnzCBayTQ2Tj2gZGIoTwNCU5g8hSSt2FadjRgeZ5Ps4LSeSx7pj+z+hImjco1TAcohwZPm84AHOGzDE4EuFJSpIg7gP+D3hLa31AKdUQkJ4yF3j3p78BePbW1pIcPFx4nXCjQxAeqNgEobXeqZR6Dgg1Pz8AvOvswARs+9d0sdjNLesbHIkw2vNdni9+JSEcrNg+CKXUQCABWGZ+Hq6U+sXZgXm6w6kXLNPSMS2EMEJJOqlfBdoDZwC01glAQyfG5PFyteb+z2IBmHxfJ2ODEW7hjh/u4I4f7jA6DOFhStIHka21PpvvV6z71Agvh6xLeTetH2RIDI0by93a7iQyJNLoEIQHKkmC2K6UuhvwVkpdDzwOrHNuWJ5rXtx+y/QXD3UzrHmpUyc5c3En0Z2ijQ5BeKCSNDE9BjTHNNzo98BZ4ElnBuXJ/rfSNErcf8d2JaxWZYOjEUJ4spJcxZQOvGR+CCf6yGp86Ya1qxgYidRicjeDZg8C4Je75PoQ4TolaWISLpCrNcsSDgHw3RM3GxwNhIfLdffupEfDHkaHIDyQJAg38c3viQBUDahIzSr+BkcjCcLdPNHxieJXEsLBStIHIVxg37GzAEz7v24GR2KSnp5Oenq60WEIIQxU1IBBH1PE5axa68edEpEH2pKUyqb9Kfj5eBNUydfocIDLY1JLH4R7uGXWLQAsHbHU4EiEJymqiUnG/nSRDxdtAWBUd6n5L+wb2Hhg8SsJ4WBFDRg0w5WBeKrftx3m+JkMalf1Z3AHuUFd2Deu3TijQxAeqKgmpiKvp9NaD3J8OJ7nvZ8TAHiwZzODIxFCCFtFNTFFAoeA2cBfgFSMc7A+byy2THe7QcZgEoXr+W1PAFbeu9LgSIQnKSpB1AF6AXcBdwOLgdla6x2uCKy8s04O7nDfg3Bvw5oPMzoE4YGK6oPIwVTie5lSyhdToohVSr2utf64uA0rpRoA32JKNLnANK31FMeEXbY9+c1ay/TTA1u5xX0Pwr2NaTvG6BCEByryRjlzYuiPKTmEAVOBBSXcdjbwjNZ6s1KqMrBJKbVCa73zKuIt87JzctmVfAaA14ZF0FHGmRZCuKmiOqlnAC2ApcBrWuvtpdmw1voocNQ8fV4ptQuoD3h0gpiyeBsAUc3rSXIQJRY1PQqA2NGxhsYhPEtRZxAjgQtAY+Bxq7LTCtBa6xJXk1NKhQE3Yurszr9sLDAWIDQ0tKSbLJO2/XuK37YkA6amJXcmpTbcy+jw0UaHIDyQ0tq5Y/8opQKB1cBbWusim6ciIiJ0fHz5vD9Pa01f80BAD/RoytBOjQyOSAhRHiilNmmtI6gY0N4AACAASURBVJyxbafWYlJK+QDzgVnFJYfyLCsn15IcgDKRHKQWk3vJyskiKyfL6DCEh3FaNVdlapP6Ctiltf6Ps/bj7lZsSWbSL1ssz2c+XjYuaZVaTO6l18xegPRBCNdyZrnvzpj6MbYppRLM817UWi8p4jXlyufLd/DzhqTLz8d2pVbVsnFJa2SkjIHsTh5s86DRIQgP5LQEobVegwfffZ1yLtOSHAZ3aMhDvW8wNqBSatJECge6k3ta3WN0CMIDyXgQTjJiSgwA9atXKnPJASAlJYWUlBSjwxBm6VnppGdJn5BwLUkQTvDXP8ct018/EmVcIFfh119/5ddffzU6DGHWb1Y/+s3qZ3QYwsPIkKMOZH0pK8CXD99kYDSiPHk44mGjQxAeSBLEVUg9n0nsjiNMW7GrwLL+bUNpEBxoQFSiPBrWQor1CdeTBHEFlmz+11Iyw57ZT/WgeqCfCyMS5d3ZTNOY5VX9qhocifAkkiBKKCc3l3d/SmDDPyfIzMqxWda5aR16tKxPpya1sSpJIoTD3DrnVkDugxCuJQmiBE6lZTL+2/Ukp16wzJsX3ZvK/j4GRiU8yeMdHjc6BOGBJEEU48yFizzw6WrSL2XTqUltnr/9Rnx9vI0OS3iYwc0GGx2C8ECSIIpw9HQ6z85cT/qlbJ4Z1IrerRsYHZLwUCnppntSggOCDY5EeBJJEIW4lJ3Dy3M2cvJsBm/d3Z6IRjWNDkl4sCE/DAGkD0K4liSIQsz+cy//pqTxf71v8MjkILWY3Mszkc8YHYLwQJIg7EhOTePHuP30aFmf2zs0NDocQ0gtJvcysMlAo0MQHkhKbeSjtebTZTvwqeDF/Tc3NTocw0gtJvdyLO0Yx9KOGR2G8DCSIPL5658TbN6fwr03NSa4iufe7Ca1mNzL8HnDGT5vuNFhCA8jTUxWtNZMX5VIveoBDIy4xuhwDNWjRw+jQxBWnu/yvNEhCA8kCcJK/L6THDhxnuhBrang7dknVw0ayCW97qTvdX2NDkF4IM/+Fsznl/iDVA/0pXuLekaHYrhDhw5x6NAho8MQZofOHuLQWfl7CNeSBGF27Ew6G/85Qd/wBh5/9gAQExNDTEyM0WEIs5E/jWTkTyONDkN4GGliMlu59TAAfW+UphXhfiZ0m2B0CMIDSYIwW7PrKM1Dq1M7KMDoUIQooOe1PY0OQXggaUvB1Lx04MR5IhvXNjoUIezaf3o/+0/vNzoM4WHkDAJYv8c0hnRkE0kQniorK4vk5GQyMzONDsWuvJvkLgZeNDgSYRQ/Pz9CQkLw8XHdMAOSIIC4xOOEBgdSv3olo0MRBklOTqZy5cqEhYW55aBPIRdDAKjsW9ngSIQRtNakpqaSnJxMw4auK//j8U1M5zOy2HrwlDQvebjMzExq1KjhlskBTIlBkoPnUkpRo0YNl5/henyC2Lj3BLlaS/OScNvkAJCZlUlmlns2fwnXMOLz6fEJYv2e41Sr5EuT+kFGhyJEoQ6ePcjBsweNDkN4GI9OEFk5uWzcd5IOjWvh5ca/Ho3Qo0cPqcfkRupXqU/9KvWduo/Y2FgGDBhQYH5CQgJLliy5om2+/fbblumkpCRatGhxxfEJ1/PoBLH1YCrpF7Ol/8GOBg0aSD0mNxJYMZDAioFkZ2e7fN9FJYji4rFOEKLs8eirmOISj+Pr482NDWWc3/zy6jB5bJKIiio4b+hQGDcO0tOhX7+Cy0ePNj1SUmDIENtlsbHF7vKNN95g1qxZNGjQgODgYNq2bUt0dDRRUVG069CO9XHrue3W2xgyZAj3338/J0+epGbNmnzzzTeEhoYyevRoBgwYwBDzvgMDA0lLSyM2NpZXX32V4OBgtm/fTtu2bfnuu+9QSrFs2TKefPJJgoODadOmTYGYLl26xMsvv0xGRgZr1qzhhRdeYNeuXRw5coSkpCSCg4Pp3bs38fHxfPLJJwAMGDCA6Oholi1bRkZGBuHh4TRv3py33nqLnJwcxowZw7p166hfvz4LFy7E39+/2GMjjOGxZxBaa9bvOU7ba4Px9fE2Ohy3I7WYXCs+Pp758+fz999/s2DBAuLj422WHz55mC8XfMkzzzzDo48+yr333svWrVsZMWIEjz/+eLHb//vvv5k8eTI7d+5k//79rF27lszMTMaMGcOiRYv4888/OXas4IBEFStW5PXXX2fYsGEkJCQwbNgwADZt2sTChQv5/vvvC93nu+++i7+/PwkJCcyaNQuAf/75h0ceeYQdO3YQFBTE/PnzS3OYhIt57BnE/uPnOXkuk5E3NTY6FLdkry3aoxT1iz8goOjlwcElOmOwtmbNGm699VbLr+mBA22HGL337nsJqWK6FyIuLo4FCxYAMHLkSJ599tlit9++fXtCQkyvDw8PJykpicDAQBo2bMj1118PwD333MO0adNKFO+gQYOu6Jd/w4YNCQ8PB6Bt27YkJSWVehvCdZx2BqGU+lopdUIptd1Z+7ga8ftOABDRqKbBkbin4OBggoOl6c1VtNZFLg8OCqZSRfs3cuZd/lihQgVyc3Mt27t06ZJlHV9fX8u0t7e3pe/gSi+drFTpcizW+wWKvFa/sDiEe3JmE9N0wG1HOdm49ySNalehRmXPHVa0KImJiSQmJhodhsfo0qULixYtIjMzk7S0NBYvXmyzPDM7k/SsdAA6derEnDlzAJg1axZdunQBICwsjE2bNgGwcOFCsrKyitxn06ZNOXDgAPv27QNg9uzZdterXLky58+fL3Q7YWFhJCQkkJuby6FDh9iwYYNlmY+PT7FxCPfltAShtf4DOOWs7V+NC5lZ7Dh0mojr5OyhMHFxccTFxRkdhsdo164dgwYNonXr1gwePJiIiAiqVq1qWX4s7ZhlwKCpU6fyzTff0KpVK2bOnMmUKVMAGDNmDKtXr6Z9+/b89ddfNr/y7fHz82PatGn079+fLl26cM019ofZ7d69Ozt37iQ8PJy5c+cWWN65c2caNmxIy5YtiY6OtunsHjt2LK1atWLEiBGlPibCeKq4U9ur2rhSYcCvWutCL35WSo0FxgKEhoa2PXjQ+TcD/bnrKG/O28ykUZG0DK3u9P2VRdOnTwdg9OjRhsbhKrt27aJZs2aGxpCWlkZgYCDp6el069aNadOmWb5s884eAnykHL0ns/c5VUpt0lpHOGN/hndSa62nAdMAIiIinJetrGz45wSBfhW4IUTunhbuY+zYsezcuZPMzExGjRpl80tcEoMwguEJwtVytSZ+30naXlsTby+PvcpXuKGiLhm9cOkCQKEd1UI4g8d9Q+49epZTaRdpd10to0MRosSSzyWTfC7Z6DCEh3HmZa6zgTigiVIqWSn1gLP2VRpxicfxUtD+ekkQouwIrRpKaNVQo8MQHsZpTUxa67ucte2rEbfnOM0bVKdqQEWjQxGixPx9pByFcD2PamI6dto89rSM/SDKmLRLaaRdSjM6DOFhPCpBrMsbe1qqt4oy5vC5wxw+d7jA/JdffpmVK1cCMHnyZNLT04vczquvvsqkSZOK3d+PP/5Is2bN6N69+5UFXELTp09HKWVT9+unn35CKcW8efMAiIqKstSmCgsL44477rCsO2/ePMul2NOnT6dmzZqEh4dbHjt37iQpKQl/f3+b+d9++22BWEpy/PKLiooiNDTU5k742267jcDAQMC2xHlsbCxKKRYtWmRZd8CAAcSay7JERUXRpEkTS4xD8hd8NIBHXcUUl3iMa2oGUk/Gni6Wx9dicjPXVLV/E9vrr79umZ48eTL33HMPAQFXf0nsV199xWeffVYgQWRnZ1OhgmO/Nlq2bMns2bMt44/MmTOH1q1bF7p+fHw8O3bsoHnz5gWWDRs2zFJVNk9SUhKNGjUiISGhyDiu9PgFBQWxdu1aunTpwpkzZzh69Gih64aEhPDWW28VqLWVZ9asWUREOOWWhiviMQniXMYltv97mqGdrjU6lDLB0+sw5d0oWJTGjRvTqVMny/p5v/zS09P54YcfbNYt7obDCxcuMHToUJKTk8nJyWHixIk0bNiQd999lwULFrB8yXKGDx/O2bNnyc3N5YYbbmD//v2WEt9HjhzhyJEjdO/eneDgYFatWsWyZct48cUXycnJITg42PIrfefOnURFRfHvv//y5JNPFqgG+/rrr7NmzRoOHDjAoEGDaN68OYsXLyYzM5MLFy4QExPDs88+y9KlS1FKMWHCBIYNG0ZsbCyvvPIKtWvXJiEhgcGDB9OyZUumTJlCRkYGP//8M40aNSrw3rt27cqff/5JVlYWFy9eZO/evZaCfvZER0fz9ttvWyrEOsLUqVMLHL/Zs2fz9ttvo7Wmf//+vPfee3ZfO3z4cObMmUOXLl1YsGABgwcPZseOHXbXbd26NVlZWaxYsYJevXo5LH5n8Zgmpg3/5I09XcfoUMoEqcXkWsuWLaNevXps2bKF7du307dvX9q0acPff/8NQExsDM2aN2Pjxo389ddfdOjQweb1jz/+OPXq1WPVqlWsWrWKkydPMmbMGObPn8+WLVv48ccfLevu3r2b5cuXs2HDBl577bUCtZJefvllIiIimDVrFh988AFgKr0yY8YMfv/9dxYsWEBCQgJbtmxh5cqVjB8/3vKrecuWLUyZMoVt27Yxc+ZM9uzZw4YNG3jwwQf5+OOP7b53pRQ9e/Zk+fLlLFy4kEGDBhV5rIYOHcrmzZvZu3dvgWVz5861aUrKyMgAYN++fTbz//zzzyKP35EjR3juuef4/fffSUhIYOPGjfz888924+nRowd//PEHOTk5zJkzx1ISvTATJkzgzTfftLtsxIgRlhjHjx9f5HZcwWPOIOISj1M90JfG9aoWv7Kw1GFq0qSJwZEYo7QlRqzXDwgIKPXr8+oYPffccwwYMICuXbsCcN1117Fr1y7i1scxYuwIyxdR3vLCrF+/nm7dutGwYUMAqle/XFKmf//++Pr64uvrS61atTh+/LilFHhhevXqZdnGmjVruOuuu/D29qZ27drcdNNNbNy4kSpVqtCuXTvq1q0LQKNGjejdu7fl/a1atarQ7Q8fPpypU6dy9uxZPvzwwyJHovP29mb8+PG888473HLLLTbL7DUx5cVSXBOTtY0bNxIVFUXNmqZ6bSNGmI79bbfdZjeeLl26MHfuXDIyMggLCyty23l/u/xJCtyvickjziAuZecQv+8kHRvXlrGnS2jo0KEMHTrU6DA8RuPGjdm0aRMtW7bkhRdesPQtdO3alaVLlxLoH8iwgcNYs2YNa9asoVu3bkVuT2tdaCnvKym5bV34r6j6bdbb9vLysjz38vIqcj/t27dn+/btpKSk0Lhx8WO0jBw5kj/++IN///232HWvRGlr1A0fPpzHHnusxP9nXnrpJd56660rCc2lPCJBJBxIJTMrR65eKoWAgACHdHaKkjly5AgBAQHcc889REdHs3nzZgC6devG5MmT6dypMyF1Q0hNTWX37t12O2ity3JHRkayevVqDhw4AMCpU44rrNytWzfmzp1LTk4OJ0+e5I8//qB9+/ZXvd133nmnxGNY+/j48NRTTzF58uSr3m8e6+PXoUMHVq9eTUpKCjk5OcyePZubbrqp0Nd27dqVF154gbvuKtntX7179+b06dNs2bLFIbE7i0c0McXtOY6fjzfhDWsYHUqZkXc6XlRnoXCcbdu2MX78eLy8vPDx8eHzzz8HTF9Ux48fp23Htpy7eI5WrVpRq1Ytu2cHY8eO5ZZbbqFu3bqsWrWKadOmMXjwYHJzc6lVqxYrVqxwSKy33347cXFxtG7dGqUU77//PnXq1GH37t1Xtd38zUXFeeCBBwq05c+dO5c1a9ZYnn/22WfUq1fP0geR5/777y/QOZ//+L3zzjt0794drTX9+vXj1ltvLTQWpRTR0dGliv+ll14qsM0RI0ZYRuoLDg62XMJsFKeW+y6tiIgInX8s3quVqzUjJsdwQ0g1Jt7Z1qHbLs+k3Ld7SUwxXTDQJNgz+4SEiceV+3a2PUdMxfnk7mlRljUMamh0CMIDlfsEEZd4DC+lpDifKNMqVpDaYcL1yn0nddye47QIrUYVf/kPJsqus5lnOZt51ugwhIcp1wniyKkLHDyZJjfHiTLvWNoxjqUdMzoM4WHKdRNTnLk4Xye5vFWUcddWkxIxwvXKd4JIPE7DWpWpU02u5xdlm4+3j9EhCA9UbpuYzqZfYsehU3JznCgzPvroI5o3b06LFi246667yMzMBODAgQNEtIug0XWNGDZsGJcuXQLg448/pkWLFvTr188yb82aNTz99NMujTs2Nvaqqv9al8QW7qXcJghTcT7k8lZRJhw+fJipU6cSHx/P9u3bLYXfAJ577jnuHns3S9YvoVq1anz11VcAfPnll2zdupUbb7yR5cuXo7XmjTfeYOLEicXuT2tNbm6uU9+TKPvKbYKISzxGcGU/rq8rxfmuhKfXYoqaHsX0hOkAZOVkETU9iu+2fgdAelY6UdOjmLt9LmC6wihqehQLdi0AICU9hajpUSxKNA0MU9LO5ezsbDIyMsjOziY9PZ169eqhteb333/n4Xsf5tpq1zJq1CibqqJZWVmkp6fj4+PDzJkz6devH9WqVbO7/aSkJJo1a8a4ceNo06YNhw4d4rfffiMyMpI2bdpw5513kpZmGrXu9ddfp127drRo0YKxY8daahPt3buXnj170rp1a9q0acO+ffsASEtLY8iQITRt2pQRI0ZY1t+0aRM33XQTbdu2pU+fPpaqr5s2baJ169ZERkby6aefluyPIlyuXCaIi1k5xO9PoWNj+yUJRPGkFpNr1a9fn+joaEJDQ6lbty5Vq1ald+/epKamEhQUhL+vPz7ePoSEhHD4sGlkuejoaDp27MjJkyfp3LkzM2bMYNy4cUXuJzExkXvvvZe///6bSpUq8eabb7Jy5Uo2b95MREQE//nPfwB49NFH2bhxI9u3bycjI4Nff/0VMJWCeOSRR9iyZQvr1q2zVG79+++/mTx5Mjt37mT//v2sXbuWrKwsHnvsMebNm8emTZu4//77eemllwC47777mDp1qqVqsHBP5bKT+u8DKVzMypHLW6+Cp9diih0da5n28faxeR7gE2DzvKpfVZvnwQHBNs/rBBb/OTx9+jQLFy7kwIEDBAUFceedd/Ldd9/Rp08f0/KM05Z18370jBw5kpEjRwLw2muv8fjjj7N06VK+/fZbGjRowIcffoiXl+1vwGuuuYaOHTsCppLgO3fupHPnzgBcunSJyMhIAFatWsX7779Peno6p06donnz5kRFRXH48GFuv/12APz8/Czbbd++vaVkeHh4OElJSQQFBbF9+3bLwDg5OTnUrVuXs2fPcubMGUvxu5EjR7J06dJij5FwvXKZIOL2HCfAtwKtw6Q435Xy9AThaitXrqRhw4aW8QcGDx7MunXrGDFiBGfOnOHI2SNUqFCBU8mnqFevns1rjxw5wsaNG3nllVdo3749cXFxvPTSS8TExBQYtSx/2e5evXoxe/Zsm3UyMzMZN24c8fHxNGjQgFdffZXMzMwSl/nOKyGutaZ58+YFzhLOnDkjZ/ZlRLlrYsrJ1azfc5x2jWri413u3p7LjB492mMK9bmD0NBQ1q9fT3p6OlprYmJiaNasGUopunfvTsLvCTSq3ogZM2YUqAA6ceJE3njjDQAyMjJQSuHl5UV6enqR++zYsSNr1661jMyWnp7Onj17LFdPBQcHk5aWxrx58wCoUqUKISEhlj6QixcvFrmPJk2acPLkSUuCyMrKYseOHQQFBVG1alVL1VVHDh0qHKvcfYPuPnyaMxcu0Umal0QZ0qFDB4YMGUKbNm1o2bIlubm5jB07FoD33nuPKZOn0LRxU1JTU3nggQcsr8sbkvTGG28ETCWwW7ZsyebNm+nbt2+R+6xZsybTp0/nrrvuolWrVnTs2JHdu3cTFBTEmDFjaNmyJbfddhvt2rWzvGbmzJlMnTqVVq1a0alTJ44dK7wDvmLFisybN4/nnnuO1q1bEx4ezrp16wD45ptveOSRR4iMjLSUtxbup9yV+562YicLNyTxwzO9qOQnNxddqbz/yJ06dTI4Etdw93LfpzJMA/5U969ezJqiPHN1ue9ydQaRqzWrdxylTaOakhyu0p49e9izZ4/RYQizkxdOcvLCSaPDEB6mXHVSbz2YSsr5TMb0ct9fgkJcieuqX2d0CMIDlasEsWJLMgEVK9BRymuIK6C1dtura7y9vI0OQRjMiO6ActPElHo+k9jtR+jVOgQ/H/nPJErHz8+P1NRUQ/4TlkRqeiqp6alGhyEMorUmNTXV5t4TVyg3ZxDz4vaTqzW3d5ChGUXphYSEkJyczMmT7tnOn1euoyQ33Ynyyc/Pz3IzoquUiwRx4Pg5Fm5MolfrEOpKaW9xBXx8fGjY0H1/XFyXY+qDkLLfwpWc2sSklOqrlEpUSu1VSj3vjH2cOJvBaz9uorK/D/ff3NQZuxDCcD7ePpIchMs57QxCKeUNfAr0ApKBjUqpX7TWO6922xcys0g+dYGN/5zgpw1J5GrNW3e3J6iSb/EvFqIMyqssOzp8tKFxCM/izCam9sBerfV+AKXUHOBWoNAEceDEeR74NBZfH298fbzx8lJczMoxPbJN/2ZeyiEzK8fymnbX1eShXjfQIDjQiW9FCGNJghBGcNqd1EqpIUBfrfWD5ucjgQ5a60fzrTcWGGt+2gLY7pSAHCcYSDE6iBKQOB1L4nQsidNxmmitKztjw848g7B3QXmBbKS1ngZMA1BKxTvrlnFHKQsxgsTpaBKnY0mcjqOUurr6REVwZid1MtDA6nkIcMSJ+xNCCOFAzkwQG4HrlVINlVIVgeHAL07cnxBCCAdyWhOT1jpbKfUosBzwBr7WWu8o5mXTnBWPA5WFGEHidDSJ07EkTsdxWoxuVe5bCCGE+yg3tZiEEEI4liQIIYQQdrlFgnBFSY5i9t9AKbVKKbVLKbVDKfWEef6rSqnDSqkE86Of1WteMMebqJTq46r3opRKUkptM8cTb55XXSm1Qin1j/nfaub5Sik11RzLVqVUG6vtjDKv/49SapQD42tidbwSlFLnlFJPusOxVEp9rZQ6oZTabjXPYcdOKdXW/LfZa37tFdUOLyTOD5RSu82x/KSUCjLPD1NKZVgd1/8WF09h79lBcTrs76xMF7j8ZY5zrjJd7OKoOOdaxZiklEowzzfkeKrCv4OM/XxqrQ19YOrA3gdcC1QEtgA3uDiGukAb83RlYA9wA/AqEG1n/RvMcfoCDc3xe7vivQBJQHC+ee8Dz5unnwfeM0/3A5ZiuielI/CXeX51YL/532rm6WpO+tseA65xh2MJdAPaANudceyADUCk+TVLgVscGGdvoIJ5+j2rOMOs18u3HbvxFPaeHRSnw/7OwA/AcPP0f4GHHRVnvuUfAi8beTwp/DvI0M+nO5xBWEpyaK0vAXklOVxGa31Ua73ZPH0e2AXUL+IltwJztNYXtdYHgL2Y3odR7+VWYIZ5egZwm9X8b7XJeiBIKVUX6AOs0Fqf0lqfBlYARY9wf2V6APu01geLid0lx1Jr/Qdwys7+r/rYmZdV0VrHadP/xm+ttnXVcWqtf9NaZ5ufrsd0X1GhiomnsPd81XEWoVR/Z/Ov25uBec6M07yfocDsorbh7ONZxHeQoZ9Pd0gQ9YFDVs+TKfrL2amUUmHAjcBf5lmPmk/hvrY6dSwsZle8Fw38ppTapExlSgBqa62PgumDBtRygzjBdO+L9X88dzuW4LhjV9887ex4Ae7H9AswT0Ol1N9KqdVKqa7meUXFU9h7dhRH/J1rAGeskqKzjmdX4LjW+h+reYYez3zfQYZ+Pt0hQZSoJIcrKKUCgfnAk1rrc8DnQCMgHDiK6VQUCo/ZFe+ls9a6DXAL8IhSqlsR6xoWp7m9eBDwo3mWOx7LopQ2LpfEq5R6CcgGZplnHQVCtdY3Ak8D3yulqrgqHjsc9Xd2Vfx3YfsjxtDjaec7qNBVC4nHocfTHRKEW5TkUEr5YPrDzNJaLwDQWh/XWudorXOB/2E6HYbCY3b6e9FaHzH/ewL4yRzTcfMpZN6p8Amj48SUwDZrrY+b43W7Y2nmqGOXjG2zj8PjNXc4DgBGmJsJMDfZpJqnN2Fqz29cTDyFveer5sC/cwqmZpMK+eY7jHnbg4G5VvEbdjztfQcVsW3XfD5L25ni6Aemu7n3Y+q4yuukau7iGBSmNrnJ+ebXtZp+ClMbKkBzbDvc9mPqbHPqewEqAZWtptdh6jv4ANuOrPfN0/2x7cjaoC93ZB3A1IlVzTxd3cHHdA5wn7sdS/J1Qjry2GEqL9ORy52A/RwYZ19MpfJr5luvJuBtnr4WOFxcPIW9ZwfF6bC/M6azT+tO6nGOitPqmK52h+NJ4d9Bhn4+HfaFcDUPTD3yezBl65cM2H8XTKdbW4EE86MfMBPYZp7/S74P/0vmeBOxuhrAme/F/IHdYn7syNs+pvbaGOAf8795HwiFadCmfeb3EWG1rfsxdRTuxeqL3EFxBgCpQFWreYYfS0xNCUeBLEy/qB5w5LEDIjCVq98HfIK5UoGD4tyLqW057/P5X/O6d5g/C1uAzcDA4uIp7D07KE6H/Z3Nn/cN5vf+I+DrqDjN86cD/5dvXUOOJ4V/Bxn6+ZRSG0IIIexyhz4IIYQQbkgShBBCCLskQQghhLBLEoQQQgi7JEEIIYSwSxKEKJOUUrFKKacPJq+UetxcYXNWvvkRSqmp5ukopVQnB+4zTCl1t719CeFKThtyVAh3pZSqoC/X+CnOOEzX7B+wnqm1jgfizU+jgDRMNy46IoYw4G7gezv7EsJl5AxCOI35l/AupdT/zDXuf1NK+ZuXWc4AlFLBSqkk8/RopdTPSqlFSqkDSqlHlVJPm4unrVdKVbfaxT1KqXVKqe1Kqfbm11cyF4nbZXfyRQAAA1RJREFUaH7NrVbb/VEptQj4zU6sT5u3s10p9aR53n8x3az1i1LqqXzrRymlfjUXVvs/4CllGj+gq1KqplJqvjmGjUqpzubXvKqUmqaU+g341nx8/lRKbTY/8s5C3gW6mrf3VN6+zNuobj4+W83Ho5XVtr82H9f9SqnHrY7HYqXUFvN7G3Z1f1XhUa727lR5yKOwB6ZfwtlAuPn5D8A95ulYzHd/AsFAknl6NKY7QCtjKntwFvPdrsBHmIqY5b3+f+bpbpjLKABvW+0jCNMdupXM203Gzl2uQFtMd6NWAgIx3Ul7o3lZEvnG3zDPjwJ+NU+/itUYCJh++XcxT4cCu6zW2wT4m58HAH7m6euB+PzbtrOvj4FXzNM3AwlW216HqZRFMKY72X0w3Rn8P6ttVc3/XuQhj8Ie0sQknO2A1jrBPL0JU9Ioziptqol/Xil1Flhknr8NaGW13mww1ftXSlVRplHWegODlFLR5nX8MH1Jg7lOvp39dQF+0lpfAFBKLcBUBvrvkrxBO3oCN6jLA3ZVUUpVNk//orXOME/7AJ8opcKBHExF4YrTBdOXPlrr35VSNZRSVc3LFmutLwIXlVIngNqYjtkkpdR7mJLMn1f4noQHkgQhnO2i1XQO4G+ezuZyE6dfEa/JtXqei+1nNn+dmLyyxndorROtFyilOgAXConxioYGLYIXEGmVCPJiIF8MTwHHgdbm12SWYNtFlW3Of6wraK33KKXaYqrr845S6jet9eslehfC40kfhDBKEqamHYAhV7iNYQBKqS7AWa31WWA58JhSlvGCbyzBdv4AblNKBSilKgG3A6X5pX0eU5NYnt+AR/OemM8Q7KkKHNWm0tgjMVU3tbe9/LGOMG83CkjRRYwboJSqB6Rrrb8DJmEaelOIEpEEIYwyCXhYKbUOU5v5lThtfv1/MVUSBXgDU9PNVmUapP6N4jaiTUM9TsdUOfQv4EutdWmalxYBt+d1UgOPAxHmjuSdmDqx7fkMGKWUWo+peSnv7GIrkG3uWH4q32tezds2ps7sUcXE1hLYoJRKwFRN9c1SvC/h4aSaqxBCCLvkDEIIIYRdkiCEEELYJQlCCCGEXZIghBBC2CUJQgghhF2SIIQQQtglCUIIIYRd/w/Ae/2PeCHiTgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.axhline(mi,label='ground truth',linestyle='--',color='red')\n",
    "for i in range(rep):\n",
    "    plt.plot(smooth_mi_list[i,:],color='steelblue')\n",
    "    plt.axvline(6000,label='switch from MINEE to MINE',linestyle='-.',color='gray')\n",
    "    for t in range(smooth_mi_list[i].shape[0]):\n",
    "        if (smooth_mi_list[0,t]>.8*mi):\n",
    "            plt.axvline(t,label='80% reached',linestyle=':',color='green')\n",
    "            break\n",
    "#plt.title(\"Plot of MI estimates against number of iteractions\")\n",
    "plt.xlim((0,smooth_mi_list[0].shape[0]))\n",
    "plt.ylim((0,mi*1.1))\n",
    "plt.xlabel(\"number of iterations\")\n",
    "plt.ylabel(\"MI estimate\")\n",
    "plt.legend()\n",
    "plt.savefig(fig_name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
